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•  Scattering  functions 

•  Hydrogen  bonding  analysis 

-  Visualization  codes 
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2.0  PROJECTS 

2.1  Ice-phobic  Materials 

The  objective  of  this  project  was  to  bring  an  understanding  to  the  adhesive  interactions 
between  coating  and  ice  and  to  establish  a  basis  of  knowledge  that  will  direct  us  towards  the 
development  of  novel  ice-phobic  coatings.  We  approached  the  problem  via  two  routes.  First,  we 
considered  investigating  the  adhesion  of  ice  on  polyurethane,  which  is  usually  the  polymer  used 
in  the  topcoat.  More  specifically  we  were  interested  in  identifying  which  chemical  groups 
promote  the  adhesion  and  what  types  of  microscopic  conformation  or  state  of  the  molecules  on 
the  surface  make  the  coating  more  vulnerable  for  the  ice  to  adhere.  We  also  wanted  to  perform 
similar  investigations  on  other  surfaces  such  as  silicon  and  teflon  in  order  to  understand  the 
general  adhesion  mechanism  of  ice  on  the  surfaces.  Secondly,  we  planned  on  approaching  the 
problem  with  an  unconventional  idea  of  understanding  how  nature  solves  it.  It  has  been  recently 
discovered  that  antifreeze  proteins,  found  in  organisms  that  can  survive  in  extremely  cold 
environments  (e.g.  Antarctic  fish),  disrupt  ice  growth  in  a  non-colligative  manner.  These  proteins 
may  directly  be  utilized  as  ice  inhibitors  or  the  knowledge  of  how  they  function  can  be  used  for 
designing  novel  coating  systems.  Thus  we  planned  on  modeling  the  binding  of  ice  to  an 
antifreeze  protein.  This  would  provide  insight  into  how  the  ice  fonnation  problem  is  solved  in 
nature  which  in  turn  will  provide  clues  as  to  how  this  problem  may  be  solved  via  synthetic 
methods.  Therefore,  how  they  function  could  be  used  to  design  biomimetic  coating  systems  with 
inherent  anti-icing  capability. 

2.1.1  Modeling  ice  crystal 

We  started  with  reviewing  the  literature  on  the  modeling  of  ice  surface  and  analyzed 
various  types  of  crystal  forms  of  ice.  It  was  feasible  to  model  ice  lh,  which  is  the  most 
energetically  favorable  form.  Besides  the  positions  of  the  oxygen  atoms  that  has  to  satisfy  the 
experimentally  determined  unit  cell  for  ice  lh,  the  positions  of  the  hydrogen  atoms  is  also 
important.  Namely,  there  is  randomness  in  the  hydrogen  bonding.  This  randomness  has  to  satisfy 
the  “ice  rules”  determined  by  Bernal  and  Fowler.  Additionally,  the  unit  cell  has  to  have  zero  net 
dipole  moment  and  minimal  net  quadrupole  moment.  I  implemented  the  results  of  a  work  by 
Hayward  and  Refiners  so  that  all  these  requirements  are  met  without  having  to  go  through  time 
consuming  Monte  Carlo  simulation  to  create  the  unit  cells.  This  implementation  is  realized  by  a 
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Fortran  code  I  developed  using  GNU  Fortran77.  The  resulting  unit  cells  can  directly  be  used  for 
any  future  modeling  studies  that  involve  ice  lh. 


Figure  1.  (Left)  Ice  lh.  (Right)  An  anti-freeze  protein. 


2.1.2  Modeling  of  polyurethane  resin 

Polyurethane  (the  main  constituent  in  topcoat  of  an  aircraft  paint)  slab  was  modeled  in 
order  to  create  the  surface  on  which  the  ice  will  adhere  and  the  adhesive  interactions  of  will  be 
studied.  The  molecule  used  in  the  modeling  is  a  dimer  of  a  product  of  the  reaction  of  m- 
Phenylene  diisocyanate  with  ethylene  glycol.  This  molecule  has  been  chosen  for  its  simplicity  in 
order  to  evaluate  the  success  of  the  force  field  and  the  simulation  protocol.  While  AMBER 
molecular  dynamics  program  package  has  shown  that  it  is  capable  of  modeling  and  simulating 
the  amorphous  polyurethane  system  at  the  realistic  density  successfully,  it  lacks  some  features 
that  enable  conveniently  modeling  of  a  two  phase  system  such  as  ice  and  polyurethane. 
Therefore,  we  intended  incorporating  LAMMPS  molecular  modeling  package  for  modeling  the 
interface.  Thereore,  codes  needed  to  be  developed  for  the  following  purposes.  First  codes  that 
make  AMBER  coordinate  and  topology  files  portable  to  LAMMPS  environment  and  vice  versa. 
Then  programs  that  analyze  output  from  LAMMPS  simulations  and  calculate  averages  and 
standard  deviations  for  several  quantities  and  plot  the  results  in  graphs  automatically. 

Next  step  was  to  model  an  isocyanurate  molecule  and  an  oligomer  of  polyester  polyol. 
The  crosslinking  reactions  of  the  two  molecules  produce  crosslinked  polyurethane;  the  type  that 
is  commonly  used  in  aircraft  topcoats.  Then  FORTRAN  code  that  creates  a  simulation  box  of 
binary  blends  on  a  cubical  lattice  was  developed.  Using  this  initial  configuration  created  by  this 
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program  we  have  performed  equilibration  simulations  of  isocyanurate  and  polyester.  Using  this 
equilibrated  structure  I  have  built  another  simulation  box  of  the  topcoat  (isocyanurate  and 
polyester)  and  the  ice  crystal.  This  initial  configuration  was  further  equilibrated  and  the  ice- 
topcoat  interface  was  obtained.  Next,  we  began  working  on  calculating  the  energy  of  adhesion. 

For  the  purpose  of  making  crosslinked  polyurethane  thus  modeling  of  a  more  realistic 
coating  system,  I  have  also  been  working  on  the  modeling  the  crosslinking  reaction.  As  reactions 
are  not  supported  by  simulation  software,  it  was  necessary  for  me  to  develop  a  programming 
code  that  performs  this  task.  We  completed  developing  the  crosslinking  code  “xlinker”.  The 
resulting  over  1300  line  long  code  will  likely  be  an  important  tool  for  creating  crosslinked 
coating  systems  towards  the  modeling  and  simulation  efforts  in  the  paints  and  coatings  group. 
The  main  function  of  the  code  is  to  determine  the  candidate  reacting  pairs  and  carry  out  the 
reactions.  In  doing  this  it  sorts  the  reacting  pairs  with  respect  to  their  closeness  thus  closer  pairs 
are  given  priority  in  the  order  of  reaction.  From  the  modeling  point  of  view,  the  reaction  of  a  pair 
of  reacting  groups  mean  that  a  set  of  changes  occur  on  several  atoms  in  the  reacting  molecules. 
These  modifications  include  change  of  atom  types,  bond  types,  angle  types  and  dihedral  angle 
types.  Additionally,  some  existing  bonds,  angles  and  dihedrals  removed  and  new  ones  created. 

An  important  feature  of  the  code  is  that  it  performs  crosslinking  reaction  of  molecules  in 
the  primary  simulation  box  with  their  images.  Therefore,  one  can  build  crosslinked  molecules 
with  infinite  length  in  any  of  the  three  dimensions.  The  modeled  system  therefore  becomes  very 
similar  to  the  real  systems.  Figure  (a)  shows  a  two  dimensional  projection  of  a  crosslinked 
molecule  that  is  infinite  in  x  and  y  directions. 

The  crosslinked  system  needs  to  be  characterized  in  order  to  know  the  micro-scale 
architecture  of  the  crosslinked  molecules.  The  code  perfonns  a  complete  analysis  and  reports  a 
wide  set  of  information  including  density  of  crosslinking  reactions,  percentage  of  reactions, 
number  of  newly  formed  chains  and  their  size,  the  members  of  newly  formed  chains,  which 
reactive  group  in  which  molecules  is  linked  with  which  other  groups  etc. 
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Figure  2.  (Left)  Repeating  group  in  polyurethane  resin  and  the  ice  polyurethane  interface. 
(Right)  Two  dimensional  view  of  cross-linked  polyurethane  molecule. 


2.2  Tunable  and  Adaptable  Adhesion 

2.2.1  Conformation  and  Dynamics  of  Self- Assembled  Monolayers 
Introduction 

Self-assembled  monolayers  (SAMs)  of  alkanethiols  provide  many  unique  features  that 
make  them  attractive  to  various  applications,  including  electrochemical  sensors1’2,  biosensors3, 
patterning  of  surfaces4,  friction  and  lubrication  control5,  protective  coatings6,  and  molecular 
electronics7.  Foremost,  they  are  relatively  easily  formed  through  solution  or  gas  phase  deposition 
and  form  a  highly  ordered  and  densely  packed  monolayer.  Thiol  groups  chemisorb  or  physisorb 
onto  surfaces  of  Au,  Ag  and  Cu;  both  experimental  and  theoretical  evidence  has  proven  that 
sulfur  atoms  coordinate  strongly  to  a  Au  surface.  ’  .  Another  important  feature  of  SAMs  is  their 
tunability.  The  surface  properties  can  easily  be  tailored  by  alternative  tail  groups  or  chain  lengths 
as  a  tuning  parameter.  For  instance,  the  selection  of  short-chain  alkanethiols  provide  a  pathway 
for  electrochemical  reactions  to  take  place  on  the  metal  surface  in  sensor  applications,  while 
long-chain  alkanethiols  enable  blocking  of  the  metal  substrate  in  applications  such  as  corrosion 
protection.  Yet,  perhaps  the  most  appealing  feature  of  SAMs  is  that  they  provide  a  molecular 
level  control  over  the  surfaces  they  are  applied  on. 

The  attractive  features  of  SAMs  have  induced  significant  attention  that  resulted  in  a 
wealth  of  studies  addressing  the  preparation,  the  characterization,  and  the  application  aspects  as 
summarized  in  comprehensive  reviews.2’10'14  Yet,  there  remain  some  essential  questions  to  be 
fully  answered,  such  as:  what  are  the  parameters  that  control  the  growth  and  the  structure,  and 
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what  are  the  driving  forces  that  render  the  monolayer  stable  or  otherwise?  Understanding  such 
questions  will  clarify  the  mechanisms  involved  in  the  growth  and  the  structure  in  SAMs,  which 
in  turn  will  help  designing  better  and  more  targeted  SAMs.11  In  addition,  a  thorough 
characterization  of  the  precise  geometry  and  the  structure  of  SAMs  is  critical  from  a 
technological  point  of  view,  particularly  in  molecular  electronics  applications  where  the 
molecular  orientation  has  direct  consequences  on  current-voltage  relationship.15'17  Evidently, 
good  control  over  the  molecule-substrate  interface  is  crucial  in  such  applications.  Controlling  the 
interface  is  not  only  a  static  matter,  the  dynamics  is  also  involved  which  is  an  overlooked  topic 
as  compared  to  other  subjects. 

SAM  structures  are  commonly  created  using  aliphatic  and  aromatic  alkanethiols. 
Aromatic  thiols  differ  from  linear  thiols  in  that  they  are  generally  more  rigid  and  they  entail 
stronger  molecule-molecule  interactions.  Additionally,  the  tunability  is  relatively  easily 
facilitated  by  aromatic  thiols  since  the  protons  around  the  phenyl  rings  can  readily  be  substituted 
with  various  chemical  groups  to  target  certain  properties.  Nevertheless,  relatively  fewer 
experimental  and  theoretical  studies  focused  on  understanding  the  system  geometry  of  SAMs  of 
aromatic  alkanethiols,  while  a  substantial  extent  of  such  effort  was  devoted  to  aliphatic 
alkanethiols.  In  one  of  such  studies  Jung  et  al.  performed  MD  simulations  of  SAMs  of  aromatic 
thiols  to  investigate  the  ordering  and  the  conformational  behavior.  Their  model  involved  the 
assumption  of  a  deprotonated  thiol  group  forming  a  covalent  bond  at  the  3 -fold  hollow  site  on  a 
Au(lll)  surface.  The  bond  fonned  between  the  sulfur  head  group  and  three  Au  atoms  were 
modeled  using  a  harmonic  potential.  Therefore,  the  sulfur  atom  was  in  effect  anchored  to  the  3- 
fold  hollow  site  and  it  was  virtually  impossible  to  observe  it  visiting  another  site  on  Au(l  1 1).  On 
the  contrary,  Kim  et  al. 19  have  adapted  a  non-bonded  potential  to  account  for  Au-biphenyldithiol 
(BPDT)  interactions.  The  resulting  binding  energy  for  the  chemisorption  was  in  reasonable 
agreement  with  experimental  measurements.  Moreover,  their  force  field  simulations  for  the  S-Au 
distance  as  well  as  the  tilt  angle  showed  a  good  agreement  with  their  complementary  Density 
Functional  Theory  (DFT)  calculations. 

Establishing  the  gold-thiol  interaction  with  a  non-bonded  potential  releases  the  restriction 
on  the  thiol  molecule  so  that  it  is  allowed  to  hop  from  one  site  to  another;  thus,  provides  a  closer 
representation  to  the  true  dynamics  of  the  organothiols.  In  fact,  the  mobility  of  the  adsorbed  thiol 
on  Au  surfaces  has  been  demonstrated  experimentally.  For  instance,  when  a  linear 
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alkanethiol  SAM  was  placed  in  a  solution  of  aromatic  alkanethiols,  the  substitution  of  the  linear 
alkanethiols  was  observed  at  the  domain  boundaries;  implying  that  the  S-Au  bond  breaks  without 
an  external  impulse.”  More  recent  evidence  to  the  mobility  was  the  observation  that  the 
stochastic  on-off  conductivity  switching,  formerly  demonstrated  in  phenylene-ethylene 
oligomers,  exists  also  with  very  simple  SAMs  of  alkanedithiols  on  Au(lll)  surface."  This 
finding  suggested  also  that  the  bond  between  the  thiol  and  Au  breaks  quite  frequently  leading  to 
on-off  mechanism. 

Our  goal  in  this  effort  was  to  characterize  the  conformational  and  the  dynamical 
properties  of  arylthiol  SAMs  and  to  contribute  to  the  understanding  of  the  influence  of  chemical 
structure  on  such  properties.  With  this  purpose  in  mind,  we  selected  six  arylthiols:  1,4- 
benzenedimethanethiol  (BDMT),  benzylmercaptan  (BM),  biphenyldithiol  (BPDT),  thiophenol 
(TP),  4-methylbenzenethiol  (TT)  and  4-methylbenzylthiol  (XT)  (see  Figure  3).  SAMs  of  all  six 
and  similar  arylthiols  with  a  single  or  two  phenyl  groups  in  standing  up  configuration  were 
formed  earlier  and  some  experimental  and  theoretical  data  were  available.  '  ’  We  followed 
a  similar  approach  to  the  work  of  Kim  et  al.19  to  account  for  Au-thiol  interactions  by 
representing  them  using  a  non-bonded  potential. 

Model 

We  used  General  Amber  Force  Field  (GAFF)  to  describe  the  intramolecular  and 
intermolecular  interactions  of  arylthiol  molecules.  GAFF  is  an  extension  to  Amber  force  field, 
which  was  designed  specifically  for  amino  acids  and  nucleic  acids,  to  cover  general  organic 
species.  The  perfonnance  of  GAFF  parameters  has  been  successfully  tested  in  terms  of 
crystalline  state  confonnations,  conformational  energies  and  dynamical  properties  for  various 
small  organic  molecules.”  ’  In  addition,  we  performed  bulk  simulation  of  BDMT  and  TP.  The 
bulk  density  was  within  3%  of  the  experimental  data.  Au-Au  interactions  were  neglected  by 
excluding  them  from  the  neighbor  list  and  Au  atoms  were  held  steady  by  setting  the  force 
imposed  by  organic  molecules  to  zero.  Keeping  the  Au  atoms  effectively  frozen  was  justified  by 
the  work  of  Kim  et  al.19  who  investigated  the  effect  of  Au  relaxation  on  conformations  of  SAMs 
of  BPDT.  They  found  that  Au  relaxation  has  a  negligible  effect  on  the  SAM  conformations.  The 
intermolecular  interactions  between  self-assembling  molecules  and  Au  atoms  were  described  by 
the  Dreiding  Force  Field  (FF)  with  optimized  parameters.  To  describe  these  non-bonded  van 
der  Walls  interactions,  we  used  an  exponential-6  form. 
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Here,  D  is  the  depth  of  the  energy  well,  R,,  is  the  distance  between  atoms,  R  is  the 
distance  at  energy  minimum  and  C  is  a  dimensionless  constant,  were  the  default  value  of  <T=  1 2 
was  used. 

Dreiding  force  field  (FF)  parameters  were  optimized  by  Jang  et  al.  "  using  DFT  targeting 
the  interactions  between  Au  atoms  and  those  in  the  molecules.  Briefly,  they  used  ethanethiol  as  a 
model  organothiol  on  a  Au32(l  11)  cluster  with  S  initially  positioned  at  FCC  and  HCP  sites  (see 
Figure  4).  The  resulting  net  binding  energy  of  24.53  kcal/mol  for  the  most  stable  FCC  site 
compares  reasonably  well  with  the  experimental  value  of  30.11  kcal/mol.  Subsequently,  the 
optimized  geometry  and  DFT  calculated  binding  energies  were  used  in  fitting  the  FF  parameters 
for  the  Au-S  interaction.  The  non-bonded  interaction  parameters  between  Au  and  other  atoms  in 
organothiols  were  obtained  by  geometric  mean  of  the  Au  parameter  in  Universal  FF34  and  the 
other  atoms  from  Dreiding  FF  .We  used  these  parameters  directly  for  our  model  (see  Table  1). 

Our  selection  of  GAFF  force  field  for  taking  into  account  the  organic  interactions  while 
using  Dreiding  FF  for  the  organic-Au  interactions  stem  mainly  from  GAFF  being  a  more 
comprehensive  force  field  than  Dreiding  FF.  Therefore,  the  various  atom  types  in  self 
assembling  molecules  can  effectively  be  represented  while  the  accurate  representation  of  the 
improved  Dreiding  FF  of  the  organic-Au  interactions  is  facilitated.  In  order  to  verify  this 
approach  we  calculated  some  confonnational  properties  of  SAMs  of  three  linear  alkanethiols: 
CioSH,  CixSH  and  C30SH  that  are  modeled  and  simulated  using  the  same  method  as  in  this 
work.35  The  tilting  angles  were  estimated  within  2°  of  the  Grazing  Incidence  X-ray  Diffraction 
(GIXD)  measurements.  While  the  tilt  direction  estimates  were  not  as  accurate  (they  differed  as 
much  as  9°),  they  showed  qualitative  agreement  with  the  experimental  findings  on  the  chain 
length  dependency.  In  addition,  we  attempted  estimating  the  melting  temperature  of  the  SAM  of 
C10SH.  For  this  purpose,  we  performed  simulations  in  which  the  SAM  was  gradually  heated.  Our 
resulting  estimate  of  Tm=  1 10  °C  agreed  reasonably  well  with  GIXD  results  of  Tm  = —  1 00  °C.1 1 
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In  order  to  construct  the  Au(l  1 1)  slab,  an  FCC  lattice  with  a  lattice  spacing  of  4.0782  A 
was  created  and  the  111  plane  was  oriented  along  the  z-axis.  The  resulting  Au(lll)  slab 
consisted  of  a  total  of  4860  atoms  arranged  in  27  x30  x  6  in  x,  y  and  z-directions. 

The  arylthiol  molecules  were  built  using  the  LEAP  module  in  AMBER  .  Point  charges 
on  atoms  were  calculated  using  the  AM1-BCC  method  which  was  parameterized  to  reproduce 
HF/6-31G*  RESP  charges.  Experimental  findings  suggest,  on  an  ideal  surface  of  Au(lll),  that 
self-assembling  organothiols  arrange  on  a  (V3>W3)-R30o  triangular  lattice  while  S  atoms  are 
positioned  on  the  FCC  (i.e.  3-fold  hollow)  site.  A  monolayer  was  formed  by  accordingly 
placing  270  organothiol  molecules.  The  conformation  of  a  self-assembling  molecule  with  respect 
to  the  substrate  is  defined  by  three  angles:  the  tilt  angle  (6),  the  azimuthal  angle  (<f>)  and  the  twist 
angle  ( <//)  (see  Figure  5).  The  orientation  of  each  self-assembling  molecule  was  initially  set  to 
standing  up  position  by  setting  <9=0,  (p  =  90  and  <//  =  0.  Afterward  they  were  randomly  rotated 
around  the  z-axis  giving  a  random  value  to  each  y/.  A  2-dimensional  view  to  a  randomly  oriented 
initial  SAM  structure  is  shown  in  Figure  6  Figure  6  (a).  The  dimensions  of  the  simulation  box 
were  set  to  be  77.86,  74.92  and  200.00  A  in  x,  y  and  z  directions,  respectively.  Setting  a  very 
large  box  dimension  in  the  z-direction  rendered  true  3 -dimensional  periodic  boundary  conditions 
as  effectively  2-dimensional. 

Simulation  Protocol 

The  MD  simulations  were  perfonned  using  LAMMPS  (Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator).40  Lennard-Jones  and  Coulombic  interactions  for  pairs  of  organic 
atoms  were  computed  using  a  switching  function  with  inner  and  outer  cutoffs  of  14  and  16  A, 
respectively,  whereas  the  cutoff  for  the  interactions  of  Au-organic  atom  pairs  was  set  to  be  14  A. 
The  long  range  Coulombic  interactions  beyond  the  cutoff  (reciprocal  sum)  were  calculated  using 
the  particle-particle  particle-mesh  Ewald  (PPPM)  solver41  with  a  precision  value  of  1.0><1 0”4.  The 
equations  of  motion  were  integrated  using  the  Verlet  algorithm  with  a  time  step  of  1  fs.  The 
SHAKE  algorithm  was  applied  to  all  hydrogen  atoms.  Temperature  in  a  canonical  (NVT) 
ensemble  was  controlled  using  the  Nose/Hoover  thermostat.43, 44 

The  system  was  minimized  using  a  conjugate  gradient  algorithm  followed  by  gradually 
raising  the  temperature  from  0  to  900  K  where  it  was  equilibrated  for  400  ps  in  the  NVT 
ensemble.  At  this  stage  of  the  equilibration,  the  sulfur  atoms  on  the  Au(l  11)  surface  were  fixed. 
Following  a  gradual  cooling  to  300  K,  the  system  was  further  equilibrated  after  removing  the 
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constraints  imposed  on  the  S  atoms  for  50  ps  and  100  ps  in  NVT  and  NVE  ensembles, 
respectively.  A  3-dimensional  view  of  the  xz-  (upper)  and  yz-  (lower)  planes  for  the  equilibrated 
SAM  of  BPDT  is  shown  in  Figure  6  (b).  Subsequent  to  the  equilibration,  the  trajectories  were 
produced  in  NVE  and  saved  for  later  analysis  for  the  first  5  ps  and  then  500  ps  with  0.001  ps  and 
0.1  ps  increments,  respectively.  While  the  most  of  the  calculations  were  made  based  on  the 
simulations  of  500  ps,  we  ran  an  additional  much  longer  simulation  of  10  ns  only  for  TT. 

Results 

To  characterize  the  confonnations  of  arylthiols  the  probability  distribution  functions  for 
the  tilt,  twist  and  azimuthal  angles  were  calculated  and  are  presented  in  Figures  5  through  7.  The 
peak  positions  for  (f>  and  i//  were  obtained  by  Gaussian  function  fits,  whereas  for  0  either 
Gaussian  or  asymmetric  double  sigmoidal  functions  were  used,  depending  on  the  shape  of  the 
distribution  function.  The  results  for  the  peak  position  as  well  as  the  full  width  at  half  maximum 
(FWHM)  were  tabulated  and  presented  in  Table  2. 

BPDT  shows  the  narrowest  distribution  of  0  with  its  peak  at  13.0°  (see  Figure  7).  A  tilt 
angle  of  13.2°  reported  by  Kim,  et  al . 19  using  DFT  is  in  excellent  agreement  with  our  model. 
Their  model  based  on  MD  simulation,  on  the  other  hand,  yielded  a  somewhat  higher  value  of 
18.4°.  BDMT  and  BM  have  broader  distributions  with  peaks  around  8°.  TP  and  TT  have  double 
peaks:  the  first  peaks  are  around  11°  and  the  second  peaks  are  at  41°  and  31°,  respectively.  In 
order  to  quantify  the  relative  weight  of  each  peak,  we  fitted  asymmetric  double  sigmoidal  and 
Gaussian  functions  for  the  fonner  and  the  latter,  respectively.  The  subsequent  integration  of 
these  functions  showed  that  the  dominant  peak  at  low  angles  comprise  70%  and  85%  of  the 
overall  distribution  for  TP  and  TT,  respectively.  Thus,  in  both  cases  only  a  small  portion  of  the 
self-assembling  molecules  tilt  at  a  higher  degree.  We  explain  the  possible  origins  of  the  second 
peak  below. 

The  experimental  data  of  the  tilting  angle  of  the  SAMs  is  commonly  obtained  from  the 
analysis  of  the  qz-dependence  of  the  grazing  incidence  X-ray  diffraction  (GIXD)  intensity.  GIXD 
measurements  of  4-methyl-4’-mercaptobiphenyl  indicated  an  upper  limit  of  19°  for  the  tilt  angle, 
which  is  substantially  smaller  than  ~30°  for  //-alkane  thiols.45  Hence,  our  simulation  results  are 
within  the  boundary  provided  by  this  experimental  measurement.  It  is  worthy  to  note  that  the 
lower  tilting  angles  associated  with  arylthiols  suggest  that  the  van  der  Waals  dimensions  for  the 
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arylthiol  molecule  are  large  enough  to  fill  the  area  available  for  each  molecule  on  the  (V3>W3)- 
R30°  lattice  structure  by  tilting  at  a  smaller  degree  than  those  of  n-alkane  thiols.11 

In  Figure  8  the  azimuthal  angle  ((f))  does  not  show  any  appreciable  deviation  from  the 
initial  configuration,  even  though  the  distributions  are  rather  dissimilar.  Resembling  the  case  of 
9,  the  BPDT  has  a  very  narrow  distribution  while  TP  has  the  broadest  distribution  of  all.  By 
looking  into  the  behavior  of  y/,  in  Figure  9  we  observe  that  there  are  two  distinct  and  alike  peaks 
as  expected  from  a  herringbone  conformation  for  all  molecules  except  TP.  The  first  and  the 
second  peaks  are  located  (on  average)  at  53°  and  127  °.  These  peaks,  on  the  other  hand,  are 
shifted  to  64°  and  115°  for  TP.  TP  has  two  additional  peaks  appearing  at  the  low  and  high  end  of 
the  distribution  scale,  specifically  at  4°  and  176°. 

The  MD  simulation  study  by  Jung,  et  al.  investigated  the  conformational  behavior  of 
BM  and  TP  that  were  anchored  on  Au(l  1 1)  each  by  a  thiolate  bond  with  (V3xV3)-7?30°  unit  cell 
arrangement.  It  is  interesting  to  note  that  while  their  model  was  constructed  in  a  fundamentally 
different  approach  than  this  work,  many  of  the  conformational  properties  agree.  For  instance, 
both  BM  and  TP  were  found  to  fonn  a  herringbone  structure,  with  the  latter  having  a  less  defined 
arrangement.  They  found  that  TP  fonned  four  i//  peaks  at  15,  65,  116  and  174°  as  compared  to 
our  findings  of  4,  64,  115  and  176°.  Moreover,  the  tilt  angle  for  BM  was  9°  versus  our  estimate 
of  8°.  However,  their  tilt  angle  for  TP  was  25°,  which  is  in  the  range  of  n-alkanethiols,  versus 
our  estimate  of  12°. 

TP  differs  from  the  other  molecules  not  only  with  its  different  number  of  peaks  in  i//  but 
also  with  the  broad  distribution  over  the  angle  scale.  Even  the  minima  regions  are  comprised  of 
appreciable  populations  of  conformations  indicating  that  TP  does  not  have  a  well-defined 
herringbone  conformation.  In  order  to  understand  the  origins  behind  the  four-peak-distribution 
we  need  to  investigate  the  two-dimensional  ordering  in  Figure  6  (c),  more  closely.  It  appears 
that,  unlike  in  SAMs  of  BPDT  (Figure  6(b))  and  the  other  molecules,  TP  has  heterogeneous 
herringbone  structure  formation  with  three-directional  orientation.  The  first  direction  along  the 
vertical  axis  has  the  same  orientation  as  in  BPDT.  The  other  two  directions  are  rotated  by  +60° 
and  -60°  with  respect  to  the  vertical  axis.  Therefore,  there  is  a  boundary  region  separating  these 
smaller  confonnationally  unalike  domains.  In  fact,  the  secondary  peak  observed  at  higher  values 
of  9  in  TP  could  be  a  result  of  the  anomalies  at  this  boundary  region. 
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Another  focus  of  our  interest  was  to  examine  the  Au(lll)  sites  (i.e.  FCC,  HCP,  Bridge 
and  Top)  that  the  thiol  head  groups  favor.  With  this  purpose  in  mind,  we  detennined  the  Au(l  1 1) 
surface  sites  at  which  the  sulfur  atom  is  positioned  at  each  time  step  of  the  simulation.  The 
criteria  used  to  detennine  the  sulfur  position  on  the  Au(lll)  surface  sites  was  as  follows.  If  a 
sulfur  atom  was  within  0.42  A  (i.e.  %  of  the  distance  between  FCC  and  HCP  sites)  of  a  site  then 
the  sulfur  was  assigned  to  that  site.  The  cases  in  which  the  sulfur  position  was  outside  of  a  0.42 
A  limit  of  any  site  were  neglected.  Such  occurrences  comprised  less  than  0.1%  of  all.  Table  3 
presents  the  fraction  of  occurrences  for  each  site  over  all  270  molecules  and  5000  steps  of  a  500 
ps-long  simulation.  It  is  clear  that  FCC  is  the  dominant  site  for  all  molecules.  The  Bridge  site  is 
rarely  (1-2%)  occupied  by  any  of  the  self-assembling  molecules.  HCP,  on  the  other  hand,  has  a 
wide  range  of  preference  from  0  to  40%.  The  most  peculiar  behaviors  of  the  site  preference  are 
observed  in  BPDT  and  TP.  BPDT  goes  almost  exclusively  to  the  FCC  site  and  it  almost  never 
visits  either  an  HCP  or  a  Bridge  site.  Interestingly,  TP  visits  FCC  and  HCP  sites  nearly  equally. 
The  occurrence  of  the  Top  site  was  inexistent. 

In  order  to  put  these  results  into  the  perspective  of  thennodynamics,  we  reviewed 
previous  theoretical  studies.  DFT  calculations  of  single  n-alkane  thiols  on  a  Au(lll)  surface 
were  perfonned  in  order  to  understand  the  energy  levels  involved  in  positioning  of  S  on  the  Au 
sites.46  For  a  single  CH3S,  the  FCC  site  was  found  to  be  the  lowest  energy  state  followed  by  HCP 
by  an  energy  cost  of  0.10  eV.  Bridge  and  Top  positions  were  local  energy  maxima  with  costs  of 
0.40  and  0.95  eV,  respectively.  Other  thiols  with  a  higher  number  of  alkane  groups  were  found  to 
follow  the  same  trend  confirming  that  FCC  site  is  the  most  stable  among  all  four  sites  on 
Au(lll).  The  above  analysis  of  site  positioning  had  shown  that  this  order  of  energy  levels 
corresponding  to  lattice  sites  is  the  major  mechanism  for  the  SAMs  of  arylthiols.  Namely,  FCC 
is  clearly  the  most  visited  site  followed  by  HCP.  Bridge  site  was  visited  dramatically  less 
frequently  than  HCP  which  represents  characteristics  of  an  energy  maximum.  A  substantially 
higher  energy  penalty  associated  with  the  Top  site  makes  its  occurrence  so  rare  that  through  the 
course  of  our  simulation  it  has  never  been  observed.  We  revisit  the  site  positioning  issue  in  the 
discussion  section. 

The  distance  between  S  and  Au  atoms  is  a  critical  parameter  in  SAMs  that  we  can  use  to 
compare  our  results  from  higher  level  theories.  Average  distances  ( dAU-s )  between  S  atoms  on  the 
Au(lll)  surface  and  its  three  nearest  Au  atoms  are  presented  in  Table  4.  The  distances  at  each 
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time  step  (every  0.1  ps)  were  first  averaged  over  the  number  of  sulfur  atoms  followed  by 
averaging  over  the  number  of  time  steps.  Standard  deviations  obtained  from  the  second 
averaging  process  were  presented  as  means  of  error.  cIau-s  for  BPDT  from  DFT  calculations  was 
reported  as  2.64  A  by  Kim,  et  a/..19  This  result  is  in  excellent  agreement  with  our  modeling 
estimate  of  2.66  A.  When  the  sulfur  atoms  were  binned  with  respect  to  their  positions  on  the 
Au(lll)  surface  sites,  the  order  of  cIau-s  was  obtained  for  all  molecules  consistently  as  FCC  < 
HCP  <  Bridge.  This  order  follows  the  intuition  due  to  the  energy  levels  of  the  sites  mentioned 
above.  Explicitly,  the  lower  the  energy  state,  the  closer  the  Au-S  pair  is.  The  variations  of  d^u-s 
among  the  molecules  are  minimal;  in  fact,  within  the  errors. 

The  fonnation  of  SAMs  is  driven  by  the  strong  interactions  between  the  head  groups  and 
the  metal  substrate  as  well  the  inter-chain  interactions  among  self-assembling  molecules.  The 
resulting  structure  comes  about  due  the  interplay  between  the  two  interactions  types.  These 
interactions  are  commonly  described  by  the  binding  energy  (Eb)  and  the  SAM  formation  energy 
(Ef).  We  define  Eb  and  Ef  in  equations  2  and  3,  respectively,  where  Evjwmu  refers  to  the  total  van 
der  Waals  energy  of  the  SAM-Au  system,  Evdw.moi  is  the  van  der  Waals  energy  of  molecule- 
molecule  interactions,  Epoh  is  the  potential  energy  of  the  SAM-Au  system,  Egas  is  the  potential 
energy  of  a  molecule  in  the  gas  phase  and  n  stands  for  the  number  of  self-assembling  molecules. 
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Eb  ranges  between  25  and  26  keal/mol  for  the  arylthiols  (Table  5).  Apparently,  it  does  not 
show  any  conclusive  variation  among  the  molecules.  Thus,  the  molecular  structure  does  not 
significantly  affect  Eb.  The  value  of  Eb  shows  a  good  agreement  with  that  of  ethanethiol  on 
Au(lll)  obtained  by  DFT  calculations,  24.53  kcal/mol.  The  experimental  results  for  various 
organothiols  measured  by  helium  beam  reflectivity  and  temperature -programmed  desorption  also 
showed  that  the  backbone  has  little  effect  on  Eb,  which  ranged  from  29.6  to  31.8  kcal/mol.  Our 
modeling  results  are  in  reasonable  agreement  with  these  experimental  measurements.  Ef  shows, 
on  the  contrary,  a  considerable  variation  with  chemical  structure.  Obviously,  the  molecular 
structure  has  a  noticeable  effect  on  the  molecule-molecule  interactions  which  lead  to  substantial 
variations  in  Ef.  The  order  of  the  magnitudes  for  Ef  is  BPDT  >  BDMT  >  XT  >  BM  >  TT  >  TP. 
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While  our  previous  analysis  of  the  sulfur  sites  (see  Table  3)  showed  the  probability  of 
finding  self-assembling  molecules  at  a  specific  site  and  a  given  time,  it  did  not  provide 
information  about  the  dynamics  of  the  site  positioning.  In  order  to  shed  some  light  on  the 
dynamics  aspect,  we  performed  a  correlation  function  analysis.  Accordingly,  the  correlation 
function  was  defined  by  the  following  equation: 

(p,(t  +  r0)pM) 

1  <P,«)> 

where  pit)  is  the  site  position  population  operator,  which  is  equal  to  1  when  a  sulfur  atom  is 
positioned  at  the  site  i  (e.g.  i  =  FCC,  HCP  or  Bridge)  at  time  t,  and  0  otherwise.  Thus,  c{t) 
describes  the  probability  that  a  sulfur  atom  originally  positioned  at  site  i  at  t=  to  stays  at  the  same 
site  at  time  t=  to+t .  The  brackets  represent  averaging  over  S  atoms  and  reference  time.  While  we 
used  the  trajectory  of  500  ps  long  in  the  calculations,  we  will  present  the  data  only  up  to  100  ps 
taking  into  account  the  statistical  error  in  time  correlation  functions.  The  resulting  time 
correlation  functions  for  BM  at  the  three  observable  sites  were  presented  in  Figure  10.  c(t)  for 
the  Bridge  site  decays  very  rapidly.  FCC  site  is  the  most  stable  of  the  three:  less  than  80%  of  c{t) 
decayed  in  100  ps.  HCP  site  has  an  intermediate  level  of  stability  with  about  50%  of  decay  in 
c(t). 

c(t)  for  each  arylthiol  were  compared  in  Figure  1 1  in  order  to  characterize  the  relative 
stabilities  of  the  S  head  groups  at  the  three  sites.  The  first  observation  to  notice  is  that  the 
relaxation  is  not  a  first  order  process.  An  initial  relaxation  occurs  around  0.1  ps  and  at  longer 
times  a  second  process  comes  into  play.  Hence,  a  second  order  exponential  decay  can  be  a  good 
model  for  fitting  c{t).  Indeed,  we  obtained  the  lifetime  of  relaxations  for  HCP  and  Bridge  sites 
from  the  second  order  exponential  fits  as  presented  in  Table  6.  We  did  not  obtain  the  lifetimes 
for  FCC  site  due  to  the  lack  of  decay  in  the  correlation  functions  for  a  plausible  extrapolation  of 
the  fitting  within  the  accessible  time  scale. 

BPDT  clearly  shows  the  highest  level  of  FCC  site  stability.  In  fact,  in  this  scale  its  second 
relaxation  process  is  imperceptible.  TP,  on  the  other  hand,  has  the  fastest  decay,  indicating  that  it 
has  the  shortest  tendency  to  stay  at  FCC  site  followed  by  TT.  XT  and  BDM  have  nearly  the  same 
intermediate  levels  of  stability  towards  the  FCC  site.  By  evaluating  the  HCP  site  stability  in 
Figure  1 1  (b)  we  observe  that  the  trend  nearly  reverses;  the  lifetimes  (see  Table  6)  show  that  the 


14 


order  of  stability  is  TP  >  BM  >  XT  >  TT  >  BDMT  >  BPDT.  There  is  a  significant  difference 
between  the  longest  and  the  shortest  lifetimes;  551  vs.  9  ps.  In  contrast  to  the  two  other  sites,  the 
relative  variations  in  decaying  of  c(t )  for  different  molecules  are  not  present  for  the  Bridge  site: 
the  lifetimes  range  within  7-8  x  10'”  ps.  In  order  to  elaborate  on  the  nature  of  the  Bridge  site,  we 
generated  a  plot  that  traces  the  trajectory  of  a  TP  head  group.  TP  was  selected  for  its  relatively 
higher  mobility  than  other  molecules  as  well  as  its  nearly  equivalent  tendency  in  being  either  at 
FCC  and  HCP  sites.  It  is  clear  in  Figure  12,  where  multiple  hopping  is  observed,  that  the  head 
group  spends  the  majority  of  time  at  either  FCC  or  HCP  locations  with  Bridge  sites  visited 
momentarily  during  the  hopping  from  one  to  the  other  minima.  It  appears  that  the  sulfur  atoms 
are  found  to  be  at  Bridge  site  because  of  the  fact  that  it  lies  on  the  path  between  two  minima  (i.e. 
FCC  and  HCP).  Therefore,  Bridge  site  exhibits  characteristics  of  a  local  energy  maximum  along 
the  pathway  between  FCC  and  HCP  sites.  Note  also  that  the  Top  position  and  its  surrounding 
area,  as  indicated  previously,  have  never  been  visited.  This  observation  is  an  indication  that  the 
Top  site  is  the  global  maximum  on  the  surface.  The  DFT  calculations  of  methanethiol  and 
ethanethiol  support  these  conclusions.46 

The  dynamics  of  SAMs  were  characterized  using  a  weight  averaged  mean  square 
displacement  (<w2>),  in  which  (n2)was  derived  as  a  function  of  time  from  the  MD  trajectories 
through  the  relationship: 

(u2}(t)=({r(t  +  t0)-r(t0))2)  (5) 

where  r(to)  and  r(7+/0)  are  the  coordinates  of  atoms  at  reference  time  to  and  after  time  t.  Brackets 
represent  averaging  over  atoms  and  reference  time.  (w2)  was  obtained  after  averaging  over  500 
ps  using  multiple  time  origins. 

Until  about  2  ps,  all  arylthiols  exhibit  similar  dynamical  behaviors  as  seen  in  Figure  13. 
This  is  the  time  scale  that  corresponds  to  the  relaxation  times  for  the  vibrational  and  the 
rotational  motions  of  small  side  groups  such  as  thiol,  methyl  and  methyl-thiol,  or  rotations  of 
phenyl  groups  around  their  equilibrium  configurations.  Since  the  molecular  structures  of  the  self¬ 
assembling  molecules  are  fairly  similar,  it  is  expected  that  such  vibrational  and  rotational 
motions  are  alike.  However,  above  the  2  ps  time  scale,  <w2>  exhibits  entirely  different  behavior. 
For  instance,  (u2)  reaches  a  plateau  for  BPDT  within  a  few  ps,  while  TP  shows  no  sign  of 
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leveling  off  within  the  scale  of  the  plot  (i.e.,  100  ps).  Apparently,  there  is  a  very  rich  and  distinct 
dynamical  behavior  at  t  >  2  ps  for  molecules  of  similar  chemical  structures. 

Discussion 

We  observed  unique  dynamical  behaviors  of  SAMs  after  analyzing  (n2> ,  particularly  in 
the  longer  time  scale  region  of  Figure  13.  For  instance,  (u2)  for  TP  continues  to  rise  at  long 
times,  whereas  in  the  case  of  BPDT  there  is  an  extended  plateau.  Understanding  the  origins  of 
this  observation  requires  a  careful  analysis  of  the  different  motions  that  contribute  to  the  overall 
dynamics  of  SAMs.  When  examining  the  relaxation  of  the  HCP  sites,  we  notice  that  BPDT 
relaxes  with  a  lifetime  of  9  ps  (see  Table  6).  It  would  be  expected  that  this  relaxation  presents 
itself  in  (w2>,  perhaps  as  a  plateau.  However,  the  HCP -relaxation  merges  together  with  the 
vibrational  and  the  rotational  contributions  that  occur  in  the  same  time  scale,  resulting  in  an 
absence  of  a  separate  plateau  at  this  short  time  scale.  On  the  other  hand,  TP  has  an  HCP- 
relaxation  of  551  ps,  which  naturally  does  not  present  itself  in  this  time  scale;  instead  of  a 
leveling  off,  (w2>  continues  to  rise.  BDMT  is  the  second  molecule  after  BPDT  that  shows  the 
shortest  HCP-relaxation  with  r=73  ps.  Note  that  the  corresponding  plateau  to  this  relaxation  is 
also  observed. 

In  order  to  elucidate  this  emerging  relationship  between  the  relaxation  of  the  sulfur  at  the 
Au(l  1 1)  sites  and  the  overall  dynamics  of  SAMs,  we  plotted  overlaid  graphs  of  (u2)(t)  and  c(t) 

in  Figure  14  using  TT  as  an  exemplar  molecule.  Note  that,  in  this  particular  case,  we  present  the 
data  extending  over  a  much  longer  time  scales  (10"  ps  <  t  <  5 x  1 0  ps).  The  time  at  which  a 
single  exponential  function  decays  to  0.368  (i.e.  e'1)  is  its  lifetime.  Therefore,  we  used  this 
number  (horizontal  dashed  line)  as  guidance  to  the  point  t  where  c{t )  has  decayed  and  reached  its 
lifetime,  approximately.  Taking  the  HCP  site  into  account,  it  can  clearly  be  seen  that  when  c(t ) 
intersects  with  the  horizontal  dashed  line  (u2^j(t)  exhibits  a  leveling  off.  The  reason  for  the 

leveling  off  is  because  prior  to  the  relaxation  time  for  HCP  site  had  been  reached,  the  HCP  site 
relaxation  was  occurring  that  was  giving  rise  to  (u2^(t).  Subsequent  to  the  completion  of  the 

substantial  relaxation,  the  HCP  site  relaxation  no  longer  affects  dynamics  notably.  At  later  times, 
(u2^(t)  shows  a  rapid  increase  due  to  the  fact  that  the  FCC  relaxation  has  not  fully  taken  place  in 

that  time  scale.  Since  we  are  bounded  by  the  length  of  the  simulation  we  can  practically  achieve, 
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we  are  unable  to  observe  the  complete  decay  of  c(t)  for  FCC.  However,  at  much  longer  time 
scales  when  FCC  exhibits  a  complete  decay,  we  would  expect  to  recognize  a  corresponding 
plateau  in  (u2Xj(t).  Identifying  the  influence  of  the  Bridge  site  decay,  on  the  other  hand,  is 

complicated  by  the  co-occurrence  of  other  modes  of  motion  (i.e.,  vibrational  and  rotational). 
Moreover,  the  relaxation  of  this  site  contributes  to  the  overall  dynamics  very  minimally,  since 
there  is  a  small  number  of  thiol  groups  positioned  at  this  site  at  a  time. 

The  above  analysis  shows  that  the  hopping  of  the  sulfur  head  group  from  one  site  to 
another  results  in  a  significant  amount  of  translational  motion  in  the  entire  molecule  and  that 
overall  the  dynamics  of  SAMs  at  the  picosecond-nanosecond  time  scale  is  driven  by  this  motion. 
On  the  other  hand,  if  the  sulfur  atom  were  to  be  fixed  at  a  particular  site  as  in  many  of  the 
previous  MD  simulations  studies,  the  motion  of  the  molecule  would  stem  mainly  from  the 
rotational  and  vibrational  sources  leading  most  likely  to  indifferent  dynamics  of  SAMs. 
Therefore,  it  is  an  important  note  of  clarity  that  the  results  for  the  unique  dynamical  behavior 
accrue  from  our  model  that  treats  the  sulfur- Au  interaction  as  non-bonded.  We  discussed 
previously  in  Sec.  I  the  experimental  evidence  supporting  the  treatment  of  sulfur- Au  interaction 
via  a  non-bonded  potential.  In  addition,  we  present  below  indirect  verification  of  this  treatment 
in  capturing  proper  structural  and  dynamical  properties  of  actual  SAMs. 

One  key  objective  of  this  work  was  to  relate  some  of  the  modeling  results  to 
experimentally  measured  properties  of  actual  SAMs.  First,  the  degree  of  order  that  we  assess  by 
the  angle  distributions  and  the  site  positioning  can  be  compared  against  the  experimental 
measurements.  Dhirani  et  al.41  studied  the  degree  of  order  in  SAMs  of 
oligo(phenylethynyl)benzenethiols  with  1,  2  and  3  phenyl  rings  using  scanning  tunneling 
microscopy  (STM).  They  found  that  the  degree  of  order  increases  with  the  number  of  phenyl 
groups.  Our  finding  that  the  arylthiol  with  two  phenyl  groups  (i.e.  BPDT)  is  more  stable  than  the 
rest  of  the  arylthiols  with  a  single  phenyl  group  is  supported  with  this  finding.  Secondly,  we 
compare  the  degree  of  integrity,  in  some  respects  the  stability,  of  SAMs  as  characterized  by  the 
surface  coverage  to  our  data.  Pugmire  et  al.  prepared  SAMs  of  all  arylthiols  modeled  in  this 
work  except  BPDT  from  solution  or  vapor  on  Au  substrate.  They  reported  the  surface  coverages 
(i.e.,  the  degree  of  integrity)  of  these  SAMs  as  measured  by  X-ray  Photoelectron  Spectroscopy 
(XPS).  The  order  of  the  surface  coverages  for  the  SAMs  from  the  highest  to  the  lowest  was 
BDMT,  XT,  BM,  TT  and  TP.  It  is  quite  interesting  to  note  that  our  results  of  a  simple  and 
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idealized  model  for  £/ followed  the  same  order  (see  Table  5)  with  the  degree  of  integrity  of  the 
actual  samples  containing  a  variety  of  defects  and  impurities.  This  probable  correlation  suggests 
that  the  SAM  formation  energy  is  an  important  parameter  that  has  implications  on  the  structure 
of  SAMs.  This  inference  was  supported  also  by  Gooding  et  al.  who  pointed  to  the  positive 
correlation  between  the  molecule-molecule  interactions  and  the  thermal  stability. 

We  may  draw  attention  also  to  the  correlation  between  the  FCC  site  stability  and  the 
surface  coverage.  When  all  the  thiol  groups  are  at  the  lowest  energy  state  of  FCC  site,  they  can 
collectively  facilitate  the  molecule-molecule  interaction  energy.  This  type  of  ordering  will  bring 
about  a  cohesive  monolayer  that  will  likely  translate  into  a  higher  surface  coverage.  On  the  other 
hand,  consider  a  single  thiol  group  visiting  a  higher  energy  site  (e.g.  HCP  and  Bridge)  with  all  its 
neighbors  are  at  the  FCC  site.  Its  new  position  will  raise  the  energy  of  the  SAM  not  only  because 
of  being  at  a  higher  energy  lattice  site,  but  also  because  it  will  inflict  the  neighboring  molecules 
to  rearrange  their  conformations  due  to  closer  contact.  As  this  type  of  visits  become  more 
frequent  the  cohesive  nature  of  the  SAM  will  diminish  leading  perhaps  to  a  lesser  degree  of 
integrity.  We  argue  that  this  scenario  not  only  links  the  FCC  site  stability  estimates  of  our 
idealized  model  to  the  surface  coverages  of  actual  SAMs;  but  also  it  is  possibly  the  mechanism 
relating  dynamics  to  structure. 

Summary 

We  predicted  the  conformational  and  the  dynamical  properties  of  SAMs  for  a  series  of 
arylthiols  on  a  Au(lll)  surface  using  MD  simulations.  The  results  of  arylthiol  configurations 
showed  a  good  agreement  with  DFT  calculations  as  well  as  experiments.  Treating  Au-thiol 
interactions  using  a  strong  non-bonded  interaction  potential  enabled  the  thiol  groups  to  visit 
other  sites  besides  FCC.  BPDT  was  found  to  form  the  most  well-defined  herringbone  structure 
with  thiol  head  groups  positioned  predominantly  at  the  energetically  favored  FCC  site.  The 
stability  of  the  head  groups  at  each  Au  site  was  characterized  by  time-dependent  correlation 
functions.  They  revealed  that  the  FCC  site  has  the  longest  relaxation  time  followed  by  the  HCP 
site,  while  the  Bridge  site  showed  characteristics  of  a  local  energy  maximum.  Moreover,  these 
site  dynamics  exhibited  a  direct  correspondence  to  the  overall  dynamical  behaviors  of  arylthiols. 
Coupled  with  our  observation  that  there  is  a  correlation  between  the  SAM  formation  energy  (Ef) 
and  the  FCC  site  stability,  we  inferred  that  Ef  is  a  crucial  parameter  that  controls  the  structural 
and  dynamical  behaviors  of  SAMs. 
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Tables 


Table  1.  Dreiding  FF  parameters  employed  for  non-bonded  interactions  between  atoms 

32 

of  organic  molecules  and  Au  atoms  from  Jang  et  al. 


Pair 

D  (kcal/mol) 

R(A) 

Au-S 

9.033 

2.682 

Au-C 

0.064 

3.561 

Au-H 

0.041 

3.082 

Table  2.  Peak  positions  and  FWHM  values  of  angle  probability  distribution  functions. 


BDMT 

BM 

BPDT 

TP 

TT 

XT 

0 

Gor  A' 

A 

A 

G 

A 

G 

G 

1st  peak 

peak  position 

8.2 

8.0 

13.0 

11.6 

10.9 

12.3 

FWHM 

9.9 

11.9 

6.1 

15.5 

12.2 

9.3 

+ 

1st  peak 

peak  position 

0.8 

-0.5 

-1.5 

1.5 

9.8 

1.5 

FWHM 

25.1 

25.3 

7.7 

26.7 

15.3 

14.1 

V 

1st  peak 

peak  position 

52.6 

53.3 

52.5 

64.  2 

53.1 

53.2 

FWHM 

18.2 

22.5 

20.4 

26.2 

23.5 

20.1 

2nd  peak 

peak  position 

127.5 

127.0 

127.5 

114.9 

127.2 

127.0 

FWHM 

18.0 

22.3 

20.3 

24.0 

24.0 

20.1 

*  G  and  A  refers  to  the  fitting  function  used.  G:  Gaussian,  A:  Asymmetric  double  sigmoid 
function. 
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Table  3.  Fraction  of  S  atoms  positioned  at  the  three  Au(l  1 1)  lattice  sites. 


FCC 

HCP 

Bridge 

BDMT 

0.93 

0.05 

0.01 

BM 

0.82 

0.16 

0.02 

BPDT 

0.98 

0.00 

0.01 

TP 

0.49 

0.40 

0.02 

TT 

0.78 

0.19 

0.02 

XT 

0.84 

0.14 

0.02 

Table  4.  Average  distance  (cIau-s)  between  S  atoms  on  the  Au(lll)  surface  and  their 
nearest  three  Au  atoms  in  A. 


Site 

BDMT 

BM 

BPDT 

TP 

TT 

XT 

FCC 

2.660  ±0.003 

2.661  ±0.004 

2.661  ±0.003 

2.667  ±  0.005 

2.663  ±  0.004 

2.662  ±0.003 

HCP 

2.709  ±0.028 

2.668  ±0.017 

2.688  ±0.081 

2.682  ±  0.006 

2.695  ±0.012 

2.689  ±0.020 

Bridge 

2.712  ±0.049 

2.721  ±0.047 

2.757  ±0.037 

2.744  ±  0.034 

2.766  ±0.029 

2.719  ±0.044 

all 

2.663  ±  0.003 

2.668  ±0.006 

2.663  ±  0.003 

2.683  ±  0.004 

2.674  ±0.004 

2.667  ±0.005 

Table  5.  Binding  (E/,)  and  formation  (Ej)  energies  in  kcal/mol. 

£(kcal/mol)BDMT  BM  BPDT  TP  TT  XT 
Eh  -25.8  -25.3  -25.1  -25.2-25.3-25.7 

Ef  -70.4  -67.1  -81.6  -61.4-64.2-68.1 
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Table  6.  Lifetime  (in  ps)  of  sulfur  atom  positions  at  the  HCP  and  Bridge  sites. 


BDMT 

BM 

BPDT 

TP 

TT 

XT 

HCPf 

73 

219 

9 

551 

190 

216 

Bridge ^ 

8  x  10'2 

7  x  10'2 

7  x  10'2 

8  x  10'2  7  x  10'2 

7  x  10'2 

T2  was 

presented 

from  a 

second  order 

exponential 

decay 

function  ( 

c(t )  =  c0  +  Axe  tlT'  +  A2e  "T2  )  over  the  data  within  1-100  ps. 

*  t  was  obtained  from  the  fitting  of  a  first  order  exponential  decay  function  ( c(t)  =  c0  +  Ae  "T ). 
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Figures 


Figure  3.  Structures  and  abbreviations  of  arylthiol  molecules  modeled  in  this  work.  1,4- 
benzenedimethanethiol  (BDMT),  benzylmercaptan  (BM),  biphenyldithiol  (BPDT),  thiophenol 
(TP),  4-methylbenzenethiol  (TT)  and  4-methylbenzylthiol  (XT). 


Figure  4.  Schematic  presentation  of  the  lattice  sites  on  Au(l  1 1)  surface.  Light  gray,  dark 
gray  and  black  beads  represent  the  Au  atoms  on  the  first,  the  second  and  the  third  layers, 
respectively. 
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with  respect  to  the  substrate. 
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Figure  6.  (a)  Views  of  xv-plane  at  (a)  initial  and  (b)  equilibrated  configurations  for  SAM 
of  BPDT  and  (c)  TP  on  Au(l  1 1)  surface.  Three  areas  are  highlighted  based  on  the  orientation  of 
the  herringbone  structure:  0°,  60°  and  -60°  orientations  are  displayed  on  white,  light  gray  and 
dark  gray  backgrounds,  respectively.  The  BPDT  atoms  that  are  used  in  producing  the  image  are 
those  enclosed  by  the  dashed  ellipse  in  Figure  3.  Views  of  (d)  xz-  and  (e)  vz-planes  for  the 
equilibrated  SAM  of  BPDT  on  Au(l  1 1). 
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Figure  7.  Probability  distribution  functions  for  the  tilt  angle  (6).  (a)  BDMT  (line),  BM 
(dashed  line)  and  BPDT  (filled  circles),  (b)  TP  (filled  triangles),  TT  (open  circles)  and  XT  (open 
triangles).  Dashed  line  in  the  inset  shows  Gaussian  function  fits  to  the  TT. 
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Figure  8.  Probability  distribution  functions  for  the  azimuthal  angle  {(/)).  (a)  BDMT  (line), 
BM  (dashed  line)  and  BPDT  (filled  circles),  (b)  TP  (filled  triangles),  TT  (open  circles)  and  XT 
(open  triangles). 
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Figure  9.  Probability  distribution  functions  of  the  twist  angle  ( y, /).  BDMT  (line),  BM 
(dashed  line),  BPDT  (fdled  circles),  TP  (filled  triangles),  TT  (open  circles)  and  XT  (open 
triangles).  The  dashed  line  with  diamond  symbols  represents  the  initial  random  distribution  of  y/. 


Figure  10.  Time  correlation  functions  (c(t))  of  sulfur  atom  positions  for  the  SAM  of  BM 
at  the  three  observable  sites  on  Au(lll).  FCC  (line),  HCP  (filled  circle)  and  Bridge  (open 
square). 
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Figure  11.  Time  correlation  functions  (c(t  j)  of  sulfur  atom  positions  for  SAMs  at  (a) 
FCC,  (b)  HCP  and  (c)  Bridge  sites  of  Au(lll)  surface.  BDMT  (line),  BM  (dashed  line),  BPDT 
(fdled  circles),  TP  (fdled  triangles),  TT  (open  circles)  and  XT  (open  triangles). 
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X  (A) 

Figure  12.  Trajectory  of  a  TP  head  group  over  500  ps  of  simulation.  FCC,  HCP,  Bridge 
and  Top  sites  were  represented  by  circles,  squares,  triangles  and  pluses,  respectively.  The  initial 
and  final  coordinates  are  indicated. 
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Figure  13.  Weight  averaged  mean  square  displacement  )  of  self-assembling 


molecules  versus  time.  BDMT  (line),  BM  (dashed  line),  BPDT  (filled  circles),  TP  (filled 
triangles),  TT  (open  circles)  and  XT  (open  triangles). 
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10'2  10'1  10°  101  102  io3  104 

t  (PS) 

Figure  14.  Overlaid  plots  of^r^(t)  (thick  line)  and  c,(t)  (thin  lines)  for  TT.  The  curves 

were  obtained  by  combining  the  plots  computed  using  trajectories  collected  at  different  intervals 
(i.e.  t  <  0.5  ps  from  0.001  ps,  0.5  <  t  <  10  from  0.1  ps  and  t  >  10  from  10  ps  intervals).  The 
vertical  dotted  lines  represent  the  points  where  the  data  were  joined  together.  Horizontal  dashed 
line  is  y  =  0.368  (i.e.,  e’1).  Intersection  points  of  c,(t)  with  the  horizontal  line  were  encircled  and 
changes  occurring  in  (u2\(t)  at  that  time  range  were  pointed  with  arrows  in  order  to  guide  the 

eye. 

2.2.2  Probing  Adhesive  Interactions  of  Self- Assembled  Monolayers  by  Simulating  Force- 
Distance  Measurements  in  Atomic  Force  Microscopy. 

Introduction 

Forces-versus-distance  measurements  using  Atomic  Force  Microscopy  (AFM)  has 
become  instrumental  in  surface  science,  colloidal  chemistry,  materials  engineering  and  molecular 
biology.  They  provide  accurate  information  on  surface  and  intermolecular  forces  that  are 
difficult  to  obtain  by  other  means.  Self-assembled  monolayers  (SAMs)  have  become  model 
systems  for  the  force-distance  measurements  since  they  easily  form  highly  ordered  and  densely 
packed  monolayers,  thus  enabling  the  measurement  of  forces  between  chemical  groups. 
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AFM  tip  as  well  as  the  substrate  can  be  coated  by  the  self-assembly  of  thiols  on  gold. 
This  approach  allows  measuring  the  forces  between  various  chemical  groups  including  methyl, 
carboxyl,  hydroxyl,  amide  and  amino.49"53 

The  experimental  measurements  of  force-distance  curves  in  SAMs  of  thiols  on  gold 
surfaces  have  yielded  conflicting  results.  Besides  the  experimental  errors,  the  inability  to  fully 
characterize  the  tip  and  the  surface  of  the  SAMs  right  prior  to  and  during  the  experiment  could 
be  a  major  contributor  to  the  inconsistent  measurements. 

In  this  effort,  we  aimed  to  investigate  the  effect  of  SAM  configuration  on  the  adhesive 
surface  forces  by  calculating  the  force-distance  curves  obtained  with  a  gold  tip  on  SAMs  of 
dithiols  at  two  different  experimentally  observed  stable  configurations:  standing-up  or  striped 
phases.  Standing-up  phase  is  formed  when  the  molecules  stand  with  their  head  groups  in  contact 
with  the  substrate  while  the  tail  groups  are  positioned  away  from  the  substrate  along  the  normal 
to  the  surface.  In  the  case  of  a  striped  phase  both  the  head  and  tail  groups  are  in  contact  with  the 
substrate.  We  further  investigated  a  possible  disordered  phase  and  a  SAM  in  standing-up 
configuration  but  when  the  tip  carried  some  residual  dithiols.  The  former  might  fonn  if  the  SAM 
was  fonned  improperly  as  well  as  if  a  properly  fonned  SAM  is  exposed  to  some  external  field 
such  as  electrical  current.  The  latter  case  is  perhaps  more  representative  of  a  standard  force- 
distance  measurement  since  the  tip  is  exposed  to  the  SAM  surface  prior  to  performing  the 
intended  force-distance  measurement.  We  investigated  the  nature  of  the  force  curves  and  how  the 
forces  are  developed  at  each  scenario. 

Model  and  Simulation  Protocol 

We  briefly  explain  the  model  and  refer  the  reader  for  further  details  to  our  earlier  work. 
We  used  General  Amber  Force  Field  (GAFF)"  to  describe  the  intramolecular  and  intermolecular 
interactions  of  arylthiol  molecules.  Au-Au  interactions  were  neglected  by  excluding  them  from 
the  neighbor  list  and  Au  atoms  were  held  steady  by  setting  the  force  imposed  by  organic 
molecules  to  zero.  The  intermolecular  interactions  between  self-assembling  molecules  and  Au 
atoms  were  described  by  the  Dreiding  Force  Field  (FF)  with  optimized  parameters.  To 
describe  these  non-bonded  van  der  Walls  interactions,  we  used  an  exponential-6  fonn. 
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Here,  D  is  the  depth  of  the  energy  well,  Ry  is  the  distance  between  atoms,  R  is  the 
distance  at  energy  minimum  and  C is  a  dimensionless  constant,  were  the  default  value  of  £"=12 
was  used.  We  used  Dreiding  force  field  (FF)  parameters  that  were  optimized  by  Jang  et  al.  for 
the  interactions  between  Au  atoms  and  those  in  the  molecules. 

Au(l  1 1)  slab  was  created  in  an  FCC  lattice  arrangement  with  a  lattice  spacing  of  4.0782 
A  after  the  111  plane  was  oriented  along  the  z-axis.  The  resulting  Au(lll)  slab  consisted  of  a 
total  of  4860  atoms  arranged  in  27  x30  x  6  in  x,  y  and  z-directions.  The  tip  was  formed  in  the 
same  lattice  structure  with  a  10  nm  diameter  (see  right  panel  in  Figure  15). 

The  arylthiol  molecules  were  built  using  the  LEAP  module  in  AMBER  .  Point  charges 
on  atoms  were  calculated  using  the  AM1-BCC  method  which  was  parameterized  to  reproduce 
HF/6-31G*  RESP  charges. 

In  standing-up  phase,  the  experimentally  suggested  arrangement  on  an  ideal  surface  of 
Au(lll)  is  the  assembly  of  organothiols  at  f A3 x A3 ) -A3 0 0  triangular  lattice  with  S  atoms  being 
positioned  on  the  FCC  (i.e.  3-fold  hollow)  site.  A  monolayer  was  formed  by  accordingly 
placing  270  organothiol  molecules  for  the  modeling  of  the  standing-up  phase.  However,  the 
molecular  arrangement  for  the  striped  phase  is  dependent  on  the  shape,  size  and  the  chemistry  of 
the  self-assembling  molecule  unlike  in  the  standing-up  phase.  Experimental  data  for  the 
arrangement  of  molecules  in  SAM  of  BPDT  in  the  striped  phase  is  not  available.  4-methyl-4 
mercaptobiphenyl  (MMB)  is  compound  that  differs  from  BPDT  only  with  a  thiol  group  instead 
of  a  methyl  group.  Therefore,  it  is  expected  that  the  molecular  arrangement  of  SAMs  of  MMB 
molecules  is  a  good  approximation  to  that  of  BPDT.  The  structure  of  MMB  SAMs  for  the  striped 
phase  was  reported  as  conforming  a  (8  x  2^3)  unit  cell  with  dimensions  of  23.08  xlO  A  and 
containing  four  molecules  per  unit  cell.45  The  head  to  head  configuration  of  molecules  results  in 
S-S  distances  of  about  2.1  A  (about  the  equilibrium  distance  of  disulfide  bond).  Unlike  MMB, 
BPDT  has  another  thiol  group  at  the  tail  position  which  enables  the  possibility  of  the  fonnation 
of  disulfide-like  bonds  between  tail  groups  as  well  as  head  and  tail  groups.  Therefore,  we  used 
Morse  bond  potential  to  describe  the  bonds  between  neighboring  sulfur  atoms. 

i  —a(r—r  )  1^ 

l-<?  v  0,\  (2) 

D  is  the  depth  of  the  potential  well,  a  is  the  stiffness  parameter  and  rQ  is  the  equilibrium  bond 
distance.  The  parameters  are  D  =  14.0  kcal/mol,  a  =  3.44  A'1  and  rQ  =  2.1  A.54  This  approach  a 


F  =  /) 

^  Morse 
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stable  structure  that  follows  the  molecular  arrangement  of  MMB  SAM  as  shown  in  Figure  15. 
The  unit  cell  dimensions  as  well  as  tilting  of  molecular  axis  with  respect  to  vertical  axis  by  a 
small  angle  are  preserved.  Since  both  thiol  groups  are  at  contact  with  the  gold  surface,  molecules 
bend  along  their  molecular  axis.  In  addition,  the  neighboring  sulfur  atoms  are  found  at  either  at 
~2.2  or  ~4.0  A  separation.  They  are  either  at  equilibrium  separation  of  a  S-S  bond  (2.1  A)  or  at 
equilibrium  separation  of  van  der  Walls  interaction  (4.0  A).  During  the  course  of  simulations  the 
S-S  separations  alter  between  these  two  distances.  Even  though  our  approach  has  yielded  a  stable 
structure  consistent  with  that  of  MMB,  it  could  be  one  of  other  reasonable  striped  phases  for 
BPDT.  Thus,  we  consider  this  structure  as  an  approximate  striped  phase  that  yields  a  plausible 
surface  coverage  for  BPDT. 

The  dimensions  of  the  simulation  box  for  all  cases  except  the  striped  phase  were  set  to  be 
77.86,  74.92  and  200.00  A  in  x,  y  and  z  directions,  respectively.  In  the  case  of  the  striped  phase,  x 
and  y  dimension  were  set  at  69.21  and  69.93  in  order  to  accommodate  the  unit  cell  dimensions 
for  the  striped  phase  (i.e.  23.08  x  1 0  A).  Setting  a  very  large  box  dimension  in  the  z-direction 
rendered  true  3-dimensional  periodic  boundary  conditions  as  effectively  2-dimensional. 

The  MD  simulations  were  performed  using  LAMMPS  (Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator).40  Lennard- Jones  and  Coulombic  interactions  for  pairs  of  organic 
atoms  were  computed  using  a  switching  function  with  inner  and  outer  cutoffs  of  14  and  16  A, 
respectively,  whereas  the  cutoff  for  the  interactions  of  Au-organic  atom  pairs  was  set  to  be  14  A. 
The  long  range  Coulombic  interactions  beyond  the  cutoff  (reciprocal  sum)  were  calculated  using 
the  particle-particle  particle-mesh  Ewald  (PPPM)  solver41  with  a  precision  value  of  1.0X10"4.  The 
equations  of  motion  were  integrated  using  the  Verlet  algorithm  with  a  time  step  of  1  fs.  The 
SHAKE  algorithm  was  applied  to  all  hydrogen  atoms.  Temperature  in  a  canonical  (NVT) 
ensemble  was  controlled  using  the  Nose/Hoover  thermostat.43, 44 

The  equilibrated  structures  of  SAMs  in  standing-up  and  striped  phases  that  were 
equilibrated  in  the  same  manner  as  previously  explained  are  shown  in  Figure  16  (a)  and  (b).55  In 
order  to  create  a  disordered  assembly,  the  SAM  at  the  standing-up  phase  was  artificially  melted 
as  follows.  The  monolayer  first  was  sandwiched  between  gold  surfaces  with  a  ~15  A  gap  by 
reducing  the  simulation  box  size  in  z  to  27  A.  The  temperature  was  raised  to  900  K  and  the 
monolayer  was  melted  over  400  ps  of  simulation  in  NVT  ensemble.  In  addition,  in  order  to 
expedite  the  melting  the  S-Au  interaction  type  was  replaced  with  C-Au  at  this  stage. 
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Subsequently,  the  temperature  was  lowered  to  300  K  and  the  box  size  and  the  interactions  were 
reset  to  the  original  values.  The  resulting  structure  (shown  in  Figure  16  (c))  will  be  discussed  in 
the  next  section. 

The  tip  was  displaced  up  and  down  in  a  quasi-equilibrium  fashion  to  simulate  its 
approach  and  withdrawal.  Subsequent  to  displacing  the  tip  along  the  z-direction  with  a  selected 
speed  over  an  incremental  distance  of  0.1  A,  the  simulation  was  run  at  the  new  position  for  the 
same  number  of  simulation  steps  needed  for  the  incremental  displacement.  At  this  stage,  forces 
applied  on  the  tip  atoms  were  collected  and  then  averaged  to  yield  a  single  data  point.  The 
collection  of  such  data  points  at  varying  displacements  produces  the  force-distance  curves.  This 
approach  has  resulted  in  reproducible  force  measurements  as  exemplified  in  the  next  section. 

Our  model  for  SAMs  of  alkane  thiols  on  Au(lll)  involves  several  assumptions  and 
approximations  besides  the  standard  approximations  in  MD  simulations.  Here  we  briefly  discuss 
their  possible  consequences.  The  representations  of  S-Au  and  molecule-molecule  interactions  are 
critical  because  the  structure  of  SAMs  is  detennined  by  the  strong  interactions  between  the  head 
groups  and  the  metal  substrate  besides  the  inter-chain  interactions  among  self-assembling 
molecules.  The  SAM  structure  comes  about  due  the  interplay  between  the  two  interactions  types, 
which  are  commonly  described  by  the  binding  energy  and  the  SAM  formation  energy.  Our 
earlier  work55  based  on  the  same  approach  has  shown  a  good  agreement  in  the  binding  energy  to 
the  results  of  Density  functional  Theory  (DFT)  and  helium  beam  reflectivity  and  temperature- 
programmed  desorption  experiments  .  Another  approximation  is  keeping  the  Au  surface  rigid. 
DFT  studies  of  Kim  et  al.19  has  shown  that  keeping  the  Au  atoms  frozen  had  no  discernible  effect 
on  the  SAM  dynamics,  which  should  be  sensitive  to  surface  characteristics.  An  additional 
assumption  was  that  the  Au  atoms  carried  no  charge.  When  thiol  molecules  are  deposited  on  the 
Au  surface,  it  is  expected  that  the  Au  atoms  on  the  surface  get  polarized  and  be  subject  to  some 
charge  distribution.  However,  since  the  S-Au  interaction  is  represented  with  a  very  strong  non- 
bonded  interaction,  any  contribution  of  the  resulting  electrostatic  interaction  to  the  overall 
potential  should  be  negligible. 

Results 

One  factor  to  consider  in  perfonning  force-distance  measurements  is  the  penetration 
depth  of  tip  into  the  sample.  It  is  a  difficult  parameter  to  control  in  AFM  experiments  because 
after  the  sudden  jump  into  direct  contact  with  the  surface,  only  force  applied  at  is  being 
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controlled.  For  instance  whether  the  tip  would  penetrate  through  a  lipid  bilayer  has  a  probability 
that  depends  on  the  force  applied  and  the  velocity  of  the  tip.56  In  our  model  since  we  control  the 
displacement  the  tip,  we  can  study  the  effect  of  the  penetration  depth.  Therefore,  we  start  our 
analysis  of  force-distance  curves  by  varying  the  level  of  approach  into  the  SAM  in  standing-up 
phase. 

When  the  tip  was  brought  to  2.4  A  from  the  surface  in  (a),  the  approach  and  withdrawal 
curves  overlap,  because  the  tip  was  not  wetted  at  this  separation.  The  well  in  the  force  curve  is 
due  to  the  strong  favorable  interactions  between  the  sulfur  atoms  on  the  surface  and  the  gold 
atoms  in  the  bottom  layer  of  the  tip.  The  wetting  occurred  with  16  thiols  being  lifted  with  the  tip 
when  the  tip  was  brought  to  the  same  z  position  as  the  surface  of  the  SAM.  The  force  curve 
shows  another  but  a  deeper  well  corresponding  to  the  second  layer  of  atoms  in  the  tip.  The 
approach  and  withdrawal  curves  are  no  longer  overlapping  due  to  the  thiol  groups  that  are 
attached  to  the  tip.  This  can  be  interpreted  as  an  inelastic  deformation  of  the  SAM.  As  the  tip  was 
penetrated  by  2.4,  4.7  and  7.1  A,  the  layers  that  are  higher  in  the  tip  become  available  to  the 
sulfur  atoms  on  the  surface.  However,  the  effects  of  the  layers  are  clearly  visible  only  until  the 
third  one,  because  the  repulsive  interactions  between  the  atoms  in  the  lower  layers  of  the  tip  and 
those  of  SAM  dominate  over  the  attractive  interactions  as  the  tip  penetrates  deeper  into  the  SAM. 
At  each  positive  penetration  the  tip  collected  16,  36,  56  and  98  dithiol  molecules  for  the 
penetration  depths  of  0.0,  2.4,  4.7  and  7.1  A,  respectively.  One  important  observation  is  that  the 
maximum  force  (Fm)  of  withdrawal  does  not  change  with  the  penetration  depth  (or  the  number  of 
dithiols  lifted  with  the  tip)  after  0  penetration.  At  this  penetration  level  the  second  layer  at  the 
bottom  of  the  tip  is  approximately  at  the  equilibrium  interaction  distance  with  the  SAM  surface. 

We  studied  the  effect  of  the  tip  speed  (from  0.5  to  0.02  A/ps)  at  the  approach  and 
withdrawal  stages  on  a  SAM  in  standing-up  phase  (Figure  20).  During  approach  all  curves 
overlap  for  the  first  two  valleys  (A  and  B).  Later  the  effect  of  speed  presents  itself.  The  hill 
between  B  and  C  diminishes  as  the  speed  is  lowered  indicating  that  SAM  molecules  can  be 
moved  away  easily.  By  end  of  equilibrium  stage  all  curves  merge.  Withdrawal  curve  does  not 
have  velocity  dependence  in  the  range  of  speeds  we  applied.  However,  it  is  worthy  of 
mentioning  that  in  actual  force-distance  measurements  the  tip  is  moved  at  a  much  slower  rate 
than  what  we  can  practically  achieve  by  MD  simulations. 
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Figure  21  presents  the  approach  and  withdrawal  curves  for  the  lower  coverage  striped 
phase.  We  do  not  observe  four  valleys  and  peaks  as  in  the  case  of  standing-up  phase.  The  depth 
of  the  first  well  is  much  shallower  than  that  of  SAM  in  standing-up  configuration  due  to  small 
density  of  molecules  on  the  surface.  Striped  phase  has  a  surface  density  of  only  57.6  A2/molecule 
compared  to  21.6  57.6  A2/molecule  in  standing-up  phase.  Following  another  shallow  well  at  3  A 
the  force  climbs  rapidly  since  the  molecules  get  squeezed  quickly  within  a  small  space  between 
the  tip  and  the  gold  surface.  At  the  withdrawal  stage  force  goes  through  a  minimum  and  then 
levels  off  while  either  head  or  tail  groups  attach  to  the  tip  get  lifted  of  the  gold  surface.  The 
leveling  occurs  because  as  one  end  of  a  molecule  id  lifted  for  more  than  a  couple  of  angstroms, 
the  attractive  interactions  quickly  drop.  Later  at  separation  of  around  10  A,  a  deeper  well  (11  nN) 
is  observed  before  force  eventually  vanishes  at  d>  15  A.  This  secondary  well  is  a  reflection  of 
other  molecules  that  are  bonded  to  those  primarily  attached  to  the  tip  being  striped  off  from  the 
surface  as  seen  in  Figure  17(b). 

We  now  investigate  the  force-distance  curve  for  perhaps  a  hypothetical  case  of  a 
disordered  assembly  of  BPDT  molecules.  Even  tough  there  has  been  no  report,  to  the  authors’ 
knowledge,  that  a  stable  disordered  assembly  has  been  formed,  it  may  be  possible  to  obtain  such 
an  assembly.  The  S-Au  bonds,  which  have  already  been  demonstrates  as  being  unstable  in 
various  experiments,  could  be  broken  collectively  by  an  external  stimulus.  '  For  instance, 
Paramanov  et  al.  have  shown  that  the  thiol-gold  interface  in  an  SAM  of  alkanethiolates  can  be 
disrupted  by  applying  a  negative  bias  to  the  tip  that  was  held  slightly  above  the  SAM77  They 
proposed  that  thiols  groups  attached  to  and  lifted  with  the  tip  after  the  S-Au  bonds  were  broken 
via  an  oxidative  cleavage  induced  by  the  negative  bias.  After  disrupting  the  S-Au  interface 
molecules  may  rearrange  their  conformations  and  get  trapped  in  a  meta-stable  disordered 
structure. 

Before  discussing  the  features  of  force-distance  curve  for  a  disordered  assembly,  we 
examine  the  structure  of  the  assembly  fonned  as  explained  earlier.  One  realizes  upon  close 
observation  of  its  structure  that  this  assembly  is  similar  to  the  striped  phase  at  gold  surface 
because  dithiols  cover  a  large  amount  of  the  surface  by  taking  a  lying  down  conformation.  60  % 
percent  of  the  molecules  with  a  direct  contact  to  the  surface  were  found  to  be  lying-down  and  the 
remaining  had  only  one  of  their  thiol  groups  on  the  surface.  The  surface  densities  of  lying-down 
and  standing-up  molecules  are  97.2  and  145.8  A2/molecule,  respectively.  The  combined  surface 
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density  in  the  disordered  assembly  is  58.3  A2/molecule,  which  is  about  the  surface  density  of  the 
SAM  in  the  striped  phase  (i.e.  57.6  A2/molecule). 

In  the  approach  stage  force  narrowly  goes  to  attractive  region,  because  the  molecules  that 
interact  with  the  tip  have  mostly  no  direct  contact  with  the  surface.  These  molecules  do  not  fonn 
bridges  between  the  tip  and  the  surface  that  give  rise  to  the  attractive  forces  as  in  ordered  SAMs. 
However,  force  clearly  goes  deep  into  the  attractive  region  at  the  withdrawal  stage.  The  features 
that  were  seen  for  the  SAM  in  the  standing -up  configuration  are  observed  noticeably  at  A’  and 
B’.  To  understand  the  changes  in  the  molecular  conformations  at  the  approach,  equilibrium  and 
withdrawal  stages  we  calculated  the  order  parameter. 

S  =  (|cos(0)|)  (2) 

where  6  is  the  angle  between  molecular  axis  and  normal  to  the  substrate  (z-axis).  Since  the 
BPDT  molecule  is  symmetric  we  do  not  distinguish  the  head  and  tail  group  positions.  We 
analyze  the  order  parameter  for  molecules  above  those  that  are  lying  down  at  the  surface  in 
Figure  23  (a).  The  assembly  appears  initially  to  be  nearly  entirely  disordered  with  a  small  bias 
towards  the  vertical  alignment  (i.e.  S  =  0.57  which  is  close  to  S  =  0.5  for  the  complete  disorder). 
As  the  tip  was  moved  toward  the  assembly  S  increased  first  and  then  decreased  signifying  that 
the  tip  oriented  the  molecules  vertically  when  at  slightly  above  the  assembly  and  disoriented 
them  when  inserted  into  the  assembly.  Meanwhile,  S  for  the  lying-down  molecules  increases 
from  approximately  zero  indicating  some  orientation  when  the  tip  gets  to  a  close  proximity 
(Figure  23  (b)).  This  increase  shows  that  one  of  the  thiol  groups  in  some  of  these  molecules  gets 
attached  to  the  tip.  During  the  equilibrium  S  remains  constant  for  both  types  of  molecules. 
However,  as  soon  as  the  tip  is  lifted  by  a  few  Angstroms  S  returns  to  its  original  value  for  the 
lying-down  molecules.  Contrarily,  S  increases  up  to  0.75  within  7  A  and  remains  above  the 
original  value  as  the  tip  is  retracted  for  the  molecules  above  those  that  are  lying-down.  This 
result  shows  that  there  is  a  tip  induced  ordering  in  the  disordered  assembly. 

In  an  experimental  force-distance  measurement,  the  tip  is  often  brought  to  contact  with 
SAM  prior  to  intended  approach  and  withdrawal  stages.  This  leads  to  the  gold  tip  to  accumulate 
some  residual  thiols  being  present  during  the  actual  measurement.  In  Figure  24,  we  present  the 
results  of  standing-up  SAM  when  the  tip  with  dithiols  carried  after  a  prior  approach  and 
withdrawal  experiment  was  used.  This  tip  had  73  BPDT  molecules  attached  to  it.  During  the 
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approach  force  never  goes  to  the  attractive  region,  rather  it  increases  gradually  until  about  3  A 
and  then  ramps  up  quickly  in  a  similar  fashion  to  the  no  residual  case  but  about  5  A  earlier. 
Apparently,  the  residual  molecules  are  forced  into  the  monolayer  pushing  the  ordered  molecules 
aside  that  screens  the  attractive  tip-SAM  interactions.  Some  of  the  residual  molecules  are  forced 
to  the  higher  layers  in  the  tip  (see  Figure  17  (d)).  During  the  equilibration  some  of  the  molecules 
in  SAM  get  flattened  while  others  manage  to  establish  attractive  interaction  with  the  tip.  For  this 
reason,  force  shows  two  attractive  valleys  in  the  withdrawal  stage  at  B’  and  A’.  An  interesting 
outcome  of  this  measurement  is  that  the  maximum  attractive  force  is  three  times  higher  when 
there  are  no  residual  molecules  than  the  case  with  residual  molecules  on  the  tip.  We  will  discuss 
the  significance  and  possible  consequences  of  this  observation  in  the  next  section. 

One  key  goal  of  this  study  was  to  understand  how  the  structure  of  SAM  affects  the 
maximum  adhesive  force  (Fm)  as  usually  measured  by  AFM.  Fm  and  the  distance  where  it  occurs 
(dm)  are  tabulated  in  Table  2  for  four  cases  studied. 
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TABLES 


Table  7.  Distance  and  depth  of  wells  in  force  curves  for  the  SAM  of  BPDT  in  standing- 
up  configuration. 


Approach  Withdrawal 


F  (nN) 

d(  A) 

A(nN) 

d(  A) 

A 

27 

3.2 

A’ 

23 

3.3 

B 

46 

1.0 

B’ 

53 

0.9 

C 

7 

-1.5 

C’ 

42 

-1.5 

D 

-21 

-3.3 

D’ 

14 

-4.0 

Table  8.  Maximum  adhesive  forces  (Fm)  and  distance  (dm)  for  different  assemblies 
studied. 


Approach  Withdrawal 


Fm  (nN) 

dm  (-A.) 

Fm  (nN) 

dm  (-A.) 

Standing-up 

46 

1.0 

53 

0.9 

Striped 

5 

3.1 

11 

11.0 

Residual 

0 

- 

17 

3.3 

Disordered 

4 

2.9 

27 

0.6 
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FIGURES 


Figure  15.  (Left  panel)  Views  through  x-y  plane  of  the  striped  phase.  Carbon  and  sulfur 
atoms  are  shown  in  blue  and  red,  respectively.  Hydrogen  atoms  were  omitted  for  clarity.  The 
rectangle  represents  the  initial  unit  cell  with  dimensions  of  23.08  xlO  A.  The  dashed  lines 
exemplify  two  different  separations  (~2.2  and  ~4.0  A)  of  neighboring  sulfur  atoms.  (Right  panel) 
Views  of  the  tip  through  x-z  and  y-z  planes  at  the  top  and  the  bottom,  respectively. 


Figure  16.  Views  through  x-z  (left)  and  y-z  (right)  planes  for  (a)  standing-up,  (b)  striped 
and  (c)  disordered  BPDT  assemblies.  Carbon,  sulfur  and  gold  atoms  are  shown  in  blue,  red  and 
gold,  respectively.  Hydrogen  atoms  were  omitted  for  clarity. 
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Figure  17.  A  view  of  the  SAM  in  standing-up  phase  with  gold  tip. 
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Figure  18.  Force-distance  curves  at  varying  penetration  depths:  (a)  -2.4,  (b)  0.0,  (c)  2.4, 
(d)  4.7  and  (e)  7.1  A.  A  negative  penetration  depth  indicates  no  penetration.  Thin  and  thick  lines 
represent  approach  and  withdrawal  stages,  respectively. 
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Figure  19.  Views  of  the  system  after  the  withdrawal  of  the  tip  for  (a)  standing-up,  (b) 
striped,  (c)  disordered  and  (d)  standing-up  (with  residual  BPDT  molecules  on  the  tip)  assemblies. 


Figure  20.  Force-distance  curves  for  the  SAM  at  the  standing-up  phase  at  four  tip 
movement  speeds:  0.5  (thin  line),  0.15  (dashed  line),  0.05  (solid  line)  and  0.02  (solid  line  with 
circles)  A/ps.  Approach,  stationary  and  withdrawal  stages  are  shown. 
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Figure  21.  Force-distance  curves  for  the  SAM  at  the  striped  phase.  Approach,  stationary 
and  withdrawal  stages  are  shown. 


Figure  22.  Comparison  of  force-distance  curves  of  the  disordered  assembly  (thin  lines)  to 
the  SAM  at  the  standing-up  phase  (thick  lines).  Approach  and  withdrawal  stages  are  shown. 
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Figure  23.  Order  parameter  for  the  disordered  assembly  calculated  considering  the  lying- 
down  molecules  (lower  panel)  and  the  rest  (upper  panel). 


Figure  24.  Comparison  of  force-distance  curves  of  the  SAM  at  the  standing-up  phase 
obtained  with  (thin  lines)  and  without  (thick  lines)  residual  BPDT  molecules  on  the  tip. 
Approach  and  withdrawal  stages  are  shown. 
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2.3  Self-Healing  Materials 

2.3.1  Coarse-Grained  Molecular  Dynamics  Simulations  of  Self-Healing  Polymeric 

Networks 

Introduction 

One  of  the  challenges  that  material  science  is  facing  is  to  design  smart  materials  that  can 
sense  the  presence  of  a  defect  and  actively  repair  the  defected  area.  Such  self-healing  materials 
would  significantly  extend  the  lifetime  and  the  usefulness  of  manufactured  items  such  as 
polymer  matrix  composites  as  replacements  of  metals,  and  cross-linked  systems  as  adhesives  or 
coatings. 

The  approaches  addressing  the  self-healing  phenomenon  in  materials  science  can  be 
categorized  in  three  classes:  healing  with  an  embedded  liquid,  thermally  activated  solid  phase 
healing  and  healing  of  projectile  puncture  59 .  The  first  approach  involves  embedding  capsules 
containing  a  reactive  healing  liquid  (e.g.,  resin)  that  rupture  and  release  the  liquid  upon 
deformation.  While  it  has  proven  to  be  successful  in  repairing  the  damaged  area  it  has  stability 
issues  while  processing  and  the  available  liquid  payload  capacity  is  often  limited  for  various 
reasons.  The  second  approach  provides  an  improved  perfonnance  after  every  cycle  of  heating, 
but  there  is  a  requirement  of  external  heat  supply  that  makes  the  process  non-autonomic.  The 
third  approach  has  gained  recent  interest  perhaps  because  it  neither  requires  an  external 
intervention  nor  has  stability  issues.  It  has  been  shown  that  the  target  film  (based  on  ionomers 
such  as  copolymers  of  ethylene  and  methacrylic  acid)  that  was  penetrated  by  a  projectile  with  a 
high  speed  is  healed  to  an  airtight  condition  almost  instantly.  Essentially,  the  penetration  of  the 
projectile  causes  localized  heating  near  the  puncture  which  then  closes  upon  itself  to  heal.  The 
mechanism  behind  such  a  self-healing  process  after  ballistic  penetration  involves  two 
transformations.  First,  the  film  which  has  been  locally  melted  due  to  heat  dissipation  re¬ 
crystallizes  and  solidifies.  Later,  the  reordering  of  the  physical  cross-links  takes  place.  These 
physical  cross-links,  according  to  Eisenberg-Hird-Moore  model  60,  are  primary  aggregates  that 
consist  of  several  ion  pairs  known  as  multiplets.  These  multiplets  bring  several  chains  together 
bound  by  ionic  interaction  creating  physical  cross-links.  Calorimetric  studies  61  have  shown  that 
the  re-crystallization  occurs  rapidly  while  the  refonnation  of  the  physical  cross-links  is  a  rather 
slow  process.  As  a  result,  the  material  fully  recovers  after  the  latter  process  is  complete. 
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Highly  cross-linked  epoxies  or  polyurethanes  are  a  class  of  materials  that  can  particularly 
benefit  from  a  self-healing  mechanism.  While  they  possess  a  variety  of  appealing  properties,  two 
major  shortcomings  of  highly  crossl inked  materials  are  their  brittleness  and  low  elongation, 
which  are  due  to  the  limited  degree  of  freedom  available  to  the  strands  between  the  cross-linking 
points.  The  strands  extend  to  their  limits  even  with  a  small  deformation  and  cause  breaking  of 
covalent  bonds.  Since  the  covalent  bonds  are  non-refonnable,  any  damage  that  is  caused  by 
impact  or  elongation  will  be  pennanent  and  cause  formation  of  micro-cracks.  For  instance,  cyclic 
mechanical  loads  result  in  accumulation  of  the  micro-cracks  which  then  leads  to  catastrophic 
failure. 

The  important  question  is:  would  a  similar  idea  to  that  behind  the  healing  of  projectile 
puncture  in  thermoplastic  ionomers  apply  to  thennosets?  It  is  known  that  the  breaking  of 
covalent  bonds  is  responsible  for  reducing  the  impact  resistance  properties.  The  inherent 
toughness  of  the  network  is  accommodated  mainly  by  non-covalent  interactions  such  as  van  der 
Waals  and  hydrogen  bond  interactions.  Their  effectiveness  is  due  to  their  capability  of  breaking 
and  displacing  under  stress  to  absorb  impact  energy.  However,  under  high  energy  impact 
conditions  such  weak  interactions  may  not  provide  a  sufficient  level  of  strength  in  the  material. 
Therefore,  in  theory,  incorporation  of  stronger,  yet  reversible  ionic  bonds  (i.e.  two  oppositely 
charged  neighboring  beads)  into  the  polymeric  matrix  could  be  an  alternative  and  effective 
approach  to  improve  the  impact  properties.  However,  one  needs  to  realize  that  a  similar 
mechanism  explained  above  for  thermoplastic  ionomers  is  an  unlikely  means  for  the  highly 
cross-linked  materials.  First,  there  will  not  be  re-crystallization  since  the  low  degree  of  freedom 
in  the  material  prevents  crystallization.  Secondly,  the  fonnation  of  multiplets  is  hindered  also 
due  to  the  cross-linked  nature  of  the  polymer. 

For  cross-linked  systems  with  embedded  ionomers  one  could  envision  an  alternative 
mechanism  in  which  the  multiplets  are  replaced  by  ionic  pairs  and  the  extent  of  repairing  action 
is  reduced  primarily  to  the  length  scale  of  strands  between  cross-links.  A  temporary  loss  of  the 
ionic  interaction  caused  by  deformation  can  be  restored  by  the  same  or  another  pair  of  ionic 
groups.  Moreover,  the  addition  of  ionic  pairs  which  are  stronger  than  van  der  Waals  and 
hydrogen  bond  interactions,  yet  still  reversible,  could  serve  to  increase  the  bulk  material  strength 
and  the  toughness.  The  soundness  of  this  probable  mechanism  can  suitably  be  studied  using 
molecular  modeling.  Past  molecular  simulation  studies  of  ionic  polymers  focused  mainly  on  their 
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charge  transport  properties  associated  with  applications  such  as  polyelectrolyte  membranes  in 
fuel  cells  62'64.  Accordingly,  these  modeling  studies  have  primarily  concerned  with  solubility, 
conductivity  and  diffusivity  rather  than  mechanical  properties. 

Computational  modeling  of  networked  polymers  has  gained  limited  attention  in  contrast 
to  their  widespread  use  in  many  applications.  Atomistic  models  of  such  systems  include  epoxy 
resins  that  were  modeled  using  a  single  step  cross-linking  procedure  65  and 
poly(dimethylsiloxane)  (PDMS)  network  where  a  dynamical  cross-linking  procedure  was  used 
66.  Besides  these  atomistic  models  Stevens  and  coworkers  modeled  polymeric  networks  using  a 
coarse-grained  representation  with  dynamical  cross-linking  procedure  67"69.  They  investigated  the 
effects  of  the  cross-linker  functionality,  the  network  connectivity  on  the  mechanical  properties 
and  the  influence  of  bond  density  on  interfacial  fracture.  Their  work  has  shown  that  a  coarse¬ 
grained  representation  enables  simulation  of  networked  systems  at  longer  time  and  length  scales 
while  capturing  the  experimentally  observed  essential  mechanical  features. 

Inspired  by  the  success  of  using  ionic  interactions  in  the  healing  of  projectile  puncture, 
here  we  seek  a  parallel  means  of  adding  self-healing  character  to  cross-linked  polymers. 
Specifically,  we  investigate  mechanical  properties  in  a  polymer  network  with  ion  pairs 
embedded  throughout  using  a  generic  coarse-grained  representation  of  the  system.  We  explore 
the  possibility  of  enhancing  toughness  by  ionic  links  and  the  potential  underlying  mechanisms 
using  molecular  dynamics  simulations.  We  have  paid  attention  also  to  the  effect  of  the  rate  of 
deformation  applied  to  the  network. 

This  article  is  organized  as  follows.  Section  2  describes  the  details  of  the  model  and  the 
simulation  protocol.  The  results  of  our  study  are  presented  in  Sec.  3.  Particularly,  the  stress-strain 
curves  of  the  standard  and  the  ionic  networks  upon  tensile  and  shear  modes  of  deformations  were 
presented.  The  possible  mechanisms  behind  the  results  were  discussed  in  Sec.  4.  The  conclusions 
of  this  work  are  given  in  Sec.  5  followed  by  the  appropriate  acknowledgments. 

Model  and  Simulation  Protocol 

The  cross-linked  system  was  represented  using  the  bead-spring  model  based  closely  on 
the  works  of  Stevens  and  coworkers  67‘69.  We  briefly  introduce  the  bonded  and  non-bonded 
potentials  used.  The  interactions  between  the  beads  were  described  using  6-12  Lennard-Jones 
(LJ)  potential  (Equation  (1))  with  a  2.5  d  cutoff. 
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U1J(r)  =  4ul 


(1) 


where  uo  is  the  LJ  energy,  d  is  the  LJ  diameter  and  r  is  the  distance  between  two  interacting 
beads.  Both  w0  and  d  were  taken  as  unity. 

The  beads  are  bonded  via  a  quartic  bond  potential  (Equation  (2)). 

-b1)(y~b2)y2  +U0,  r<rc j 
r>rc\ 

where  y  =  r  -  A r  that  shifts  the  quartic  center  from  the  origin  and  the  parameters  as  adopted 
directly  from  the  study  by  Stevens  et  al.  67  are  k  =  1434.3  uo/d4,  b\  =  -0.759  d,  bi  =  0.0,  A r  =  rc  = 
1.5  d  and  Uo  =  67.223  u0.  The  bond  potential  is  smoothly  cutoff  at  rc.  At  distances  larger  than  rc 
the  bond  potential  is  turned  off  that  leaves  only  LJ  term  remaining.  The  bond  breakage  is 
irreversible.  The  bond  potential  parameters  result  in  maximum  bonding  force  of  156.7  uo/d 
compared  to  maximum  LJ  force  of  2.4  uo/d. 

Coulombic  interactions  for  pairs  of  ionic  beads  were  computed  with  a  cutoff  of  also  2.5 
d. 


(3) 


where  q,  is  the  charge  on  bead  i,  so  is  dielectric  constant.  As  this  is  a  coarse  grained  model  with 
each  bead  representing  many  atoms,  the  size  of  a  bead  is  much  larger  than  atoms.  Therefore, 
when  the  charges  are  placed  at  the  center  of  large  beads,  ionic  bonds  (i.e.  two  oppositely  charged 
neighboring  beads)  are  separated  by  large  distances.  Thus,  the  Coulombic  force  between  two 
neighboring  beads  falls  short  of  a  realistic  ionic  interaction.  In  order  to  bring  the  ionic  interaction 
force  to  a  more  realistic  quantity  we  set  the  dielectric  constant  (so)  to  0.02.  This  dielectric 
constant  was  selected  in  order  to  situate  the  ratio  of  the  maximum  quartic  bonding  force  to  ionic 
interaction  force  at  about  50.  This  ratio  is  simply  an  estimate  value  that  relates  the  strength  of  an 
typical  covalent  bond  to  an  ionic  bond.  In  fact,  at  r  =  1.12  d  the  ionic  force  is  3.17  uo/d  versus 
156.7  uo/d  (i.e.  ratio  of  49.4). 

The  MD  simulations  were  performed  using  LAMMPS  (Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator)  70.  The  equations  of  motion  were  integrated  using  the  Verlet 
algorithm  with  a  time  step  of  0.05  x,  where  x  is  the  LJ  time  unit.  The  long-range  Coulombic 
interactions  beyond  the  cutoff  (reciprocal  sum)  were  calculated  using  the  particle-particle 
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particle-mesh  Ewald  (PPPM)  solver  with  a  precision  value  of  1.0X10'\  Both  temperature  and 
pressure  were  controlled  using  Nose-Hoover  formalism  with  a  damping  constant  of  0. 1  x. 

The  initial  simulation  box  consisted  of  small  chains  with  two  and  three  beads  that  were 
sandwiched  between  two  rigid  layers  of  beads.  The  chains  with  three  beads  contained  one  bead 
designated  as  the  cross-linker  bead  for  subsequent  cross-linking  process.  The  layers  were  made 
of  (111)  plane  of  face-centered  cubic  (FCC)  lattice  structure  with  lattice  spacing  of  1.2.  They 
interact  with  all  beads  except  themselves.  The  system  consisted  of  a  total  number  of  103,344 
beads  of  chain,  cross-linker  and  layer  type,  where  the  ratio  of  the  two-  and  three-bead  chains  was 
determined  by  stoichiometry.  The  liquid  mixture  was  first  equilibrated  in  canonical  ensemble 
(NVT)  at  1.0  reduced  temperature  unit  ( T *),  followed  by  the  density  equilibration  in  NPT 
(isobaric-isothermal)  ensemble.  The  pressure  was  imposed  anisotropically  in  z-axis  while  no 
pressure  scaling  was  employed  in  x  and  y  directions. 

The  liquid  mixture  was  dynamically  cross-linked  after  the  equilibration.  Briefly,  a 
varying  distance  criterion  was  used  first  to  detennine  the  candidate  pairs  (i.e.  a  bead  and  a  cross- 
linker  bead  with  a  functionality  of  four)  which  then  were  sequentially  cross-linked.  The  selected 
separation  distance  ranged  from  1.1  d  at  the  first  step  until  1.35  d  at  the  final  cross-linking  step. 
Each  step  consisted  of:  1.  detennination  of  pairs,  2.  creation  of  bonds  and  3.  equilibration  of  the 
system.  The  cross-linking  (first  to  beads  in  the  layers  and  then  in  the  bulk)  was  carried  on  for  a 
total  of  40  and  60  steps,  respectively.  The  number  of  steps  for  the  first  cross-linking  stage  was 
adjusted  to  create  adequate  cross-linking  density  on  the  layers  so  that  cohesive  failure  was 
promoted.  The  second  cross-linking  stage  was  continued  until  a  sufficient  degree  of  cross-linking 
was  achieved  in  the  bulk.  Near  the  end  of  each  cross-linking  stage,  only  a  few  or  no  additional 
cross-linking  occurred.  As  a  result,  about  60,000  cross-linking  bonds  were  created  with  nearly 
5,000  of  them  formed  to  the  beads  in  the  rigid  layers,  leading  to  98%  degree  of  cross-linking. 
After  equilibration  of  the  network  at  1.0  T*,  the  harmonic  bond  potential  used  hitherto  was 
replaced  with  the  quartic  potential.  Subsequently,  the  system  was  slowly  cooled  down  (at  a  rate 
of  1.4x1  O'5  T *  per  MD  step)  to  0.3  T*,  which  is  well  below  the  estimated  glass  transition 
temperature. 

Ionic  bonds  were  incorporated  into  the  fully  cross-linked  network  subsequent  to  its 
equilibration.  The  bond  between  a  pair  of  randomly  selected  consecutive  beads  was  removed  to 
create  ionic  bonds.  In  the  case  of  an  ionic  network,  these  beads  were  assigned  opposite  charges 
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to  create  ionic  bonds  while  for  a  standard  network  they  were  left  unchanged.  Using  this 
approach,  we  created  a  new  network  with  the  same  network  density  (82%)  as  non-ionic  one  but 
with  ionic  bonds  incorporated  within.  The  fact  that  we  used  pairs  of  beads  bonded  to  each  other 
to  create  ionic  pairs  is  actually  expected  to  reflect  a  molecular  arrangement  close  to  an  actual 
system.  Namely,  one  would  anticipate  in  the  actual  synthesis  of  a  polymeric  network  that  the 
monomers  with  opposite  charges  would  be  at  close  proximity  as  the  cross-linking  reactions  take 
place. 

Two  forms  of  deformation  (tensile  and  shear)  were  applied  in  a  quasi-equilibrium 
fashion.  In  tensile  mode,  the  upper  layer  was  displaced  along  the  z-direction  with  a  selected 
speed  over  an  increment  of  0.5  d.  Subsequently,  the  simulation  was  run  at  the  new  position  for 
the  same  number  of  simulation  steps  needed  for  the  displacement.  The  data  for  the  calculation  of 
stress-strain  curves  were  collected  at  this  stage.  The  effective  deformation  rate,  therefore,  is  half 
of  the  selected  displacement  speed.  The  two  selected  speeds  0.5x10'  and  0.5x10'“  d/x  resulted  in 
strain  rates  of  0.6xl0'3  and  0.6X10’4  1/x,  respectively.  In  the  case  of  shear  deformation  the  upper 
layer  was  displaced  in  the  x-direction. 

Results 

We  begin  our  analyses  with  a  careful  look  into  the  microstructure  of  the  ionic  network 
using  a  set  of  radial  distribution  functions  (g(r))  shown  in  Figure  25.  g {f)  chain-chain  exhibits  an 
initial  sharp  peak  located  at  0.95  and  a  broader  peak  at  1.07  d  corresponding  to  the  bonded  and 
the  non-bonded  interactions  among  the  beads  of  chain  type,  respectively.  In  other  words,  the 
quartic  bond  potential  with  the  parameters  used  in  this  study  shows  a  minimum  in  the  potential 
energy  at  0.93  d  that  agrees  with  the  position  of  the  initial  peak,  while  the  energy  minimum  in 
Lennard- Jones  potential  for  the  non-bonded  interactions  is  at  1.12  d  in  line  with  the  second  peak. 
The  remainder  of  the  g if) chain-chain  has  shallower  and  dying  out  peaks  that  is  reminiscent  of 
liquids,  g (r)xiinker_xiinker  is  very  similar  to  g(r)chain.chain  suggesting  that  the  cross-linker  beads  are 
uniformly  distributed  in  the  network.  However,  the  absence  of  a  sharp  peak  near  0.93  d  signifies 
that  cross-linkers  are  not  bonded  to  one  another.  The  sharp  spike-like  peaks  in  g {f)xiinker-iayer 
indicate  a  layer  by  layer  positioning  of  cross-linker  beads  near  the  layers.  The  first  peak  is  a 
result  of  cross-linkers  directly  bonded  to  the  layers,  while  the  second  (~1.6  d)  and  third  (-2.0  d) 
peaks  are  due  to  the  cross-linkers  that  are  separated  from  layers  by  one  and  two  layers  of  beads, 
respectively.  The  intensity  of  these  peaks  diminishes  rapidly  with  r.  This  type  of  arrangement 
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shows  that  the  network  has  a  fairly  ordered  structure  near  the  layers  (r  <  ~5  d)  that  becomes 
amorphous-like  away  from  the  layers.  The  relative  positioning  of  ionic  pairs  are  presented  by 
g(r)ion-ion.  There  is  a  single  and  dominant  peak  at  0.97  d  pointing  to  the  fact  that  they  were 
created  from  pairs  of  beads  which  were  initially  bonded  (i.e.  pairs  of  cross-linker  and  chain). 

After  this  brief  examination  of  the  microstructure  we  will  now  investigate  the  material 
responses  to  different  fonns  of  deformation  for  both  control  and  ionic  networks  starting  with  the 
tensile  deformation  in  Figure  26  (lower  panel).  There  is  a  clear  peak  at  0.1  seen  for  both 
networks  corresponding  to  the  yield  behavior.  The  yield  stress  originates  from  the  LJ 
interactions.  For  this  reason,  the  yield  peak  appears  to  be  very  distinct  due  to  the  fact  that  the 
networks  have  a  lower  degree  of  cross-linking  than  a  fully  cross-linked  system.  The  yield  point, 
which  is  effectively  independent  of  the  degree  of  cross-linking,  becomes  larger  than  the  rest  of 
the  stress-strain  curve  since  the  degree  of  cross-linking  was  lowered  to  82%.  Note  also  that  the 
equilibrium  distance  of  the  LJ  potential  is  at  1.12  d,  while  the  maximum  force  occurs  at  1.25  d. 
Therefore,  it  takes  a  strain  of  about  0.1  to  bring  the  network  from  the  equilibrium  LJ  separation 
to  a  separation  of  the  maximum  resistance. 

In  the  standard  network,  subsequent  to  the  yield  stress  there  is  plateau  region  that  lasts 
until  about  e  =  1.0.  During  this  period  the  strands  connecting  two  cross-linkers  are  being 
extended  while  bonds  remain  mostly  intact  as  described  by  Stevens  .  The  bonds  start  breaking 
appreciably  only  after  about  s  =  1.2  (see  upper  panel  in  Figure  26)  and  the  network  fails 
ultimately  at  about  s  =  2.0.  The  ionic  network  exhibits  a  similar  behavior  to  the  standard 
network,  but  with  some  imperative  differences.  First,  after  the  yield  behavior  the  values  of  stress 
do  not  drop  as  extensively  as  the  standard  network.  The  stress-strain  curve  remains  above  that  of 
the  standard  network  indicating  that  the  addition  of  ionic  bonds  improves  the  toughness  (i.e.  the 
area  under  this  curve).  Another  key  result  is  that  the  ionic  network  requires  nearly  twice  as  many 
bonds  to  be  broken  to  fail  compared  to  the  standard  network.  We  also  observe  that  the  ultimate 
failure  does  not  occur  abruptly  due  to  lower  degree  of  cross-linking  and  limited  system  size. 

Snapshots  of  the  network  under  tensile  deformation  at  five  consecutive  strains  are 
presented  in  Figure  27  to  aid  visualizing  the  defonnation  in  tensile  mode.  At  s  =  0.5,  there  are 
voids  formed  in  the  network.  They  actually  emerge  as  strain  passes  through  the  yield  point.  Their 
formation  does  not  require  breaking  of  bonds  since  they  are  already  present  but  in  compacted 
form.  Moreover,  their  presence  indicates  that  the  cross-linking  is  not  entirely  homogenous:  there 
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are  regions  of  highly  dense  degree  of  cross-linking  and  those  of  lesser  amount  of  cross-links. 
They  basically  act  as  defects  that  initiate  cracks.  At  s  =  1.0,  the  voids  grow  further  and  the 
surrounding  areas  approach  their  limits  of  the  stress  they  can  carry  without  breaking  of  bonds. 
Later,  at  e  =  1.5,  the  stretching  strands  that  bear  the  load  become  more  visible.  At  this  stage, 
there  is  a  considerable  degree  of  bond  breaking  and  the  network  starts  to  break  into  two  pieces. 
Finally,  at  s  =  2.0,  the  network  has  nearly  failed  with  only  a  single  strand  holding  the  newly 
created  interfaces. 

The  response  to  defonnation  in  shear  mode  is  presented  in  Figure  28.  The  initial  rapid 
rise  of  stress  and  the  subsequent  slowing  down  at  s  =  0.1  is  due  to  the  yield  behavior.  The 
magnitude  of  the  yield  stress  for  the  standard  network  is  0.07  uo/d3  which  is  in  agreement  with 
0.08  ug/d3  found  in  the  earlier  study  67 .  The  yield  stress  in  this  deformation  mode  stems  from  the 
frictional  resistance  between  two  sliding  surfaces  interacting  via  LJ  beads.  After  the  yield  point, 
stress  continues  to  rise  first,  at  a  lower  rate,  and  later  (s  >  -1.2)  at  a  higher  rate  until  it  reaches 
the  maximum  point  (i.e.  where  the  ultimate  failure  occurs).  Note  that  the  bond  breaking 
continues  even  after  the  maximum  point,  because  the  failure  results  in  a  rough  interface.  This 
interface  is  smoothened  while  being  sheared,  thus  bond  breaking  continues  until  two  sufficiently 
smooth  surfaces  are  attained.  Subsequently,  the  shearing  would  constitute  shearing  of  two  flat 
surfaces  interacting  via  LJ  interactions  similar  to  the  case  at  the  yield  point.  Therefore,  if 
simulation  was  pursued  further  the  stress  would  eventually  level  off  at  about  the  yield  stress. 

In  the  earlier  studies  ’  the  maximum  point  was  followed  by  a  sudden  drop  in  stress. 
This  drop  is  more  gradual  in  our  study  possibly  because  of  the  lower  degree  of  cross-linking 
since  densely  cross-linked  networks  fail  more  abruptly.  The  comparison  of  networks  shows  that 
stress  increases  in  ionic  network  more  rapidly  than  the  standard  network  and  maintains  its 
elevated  level  until  the  point  of  maximum  stress.  Similar  to  what  was  found  in  the  defonnation  in 
tensile  mode,  the  incorporation  of  ionic  pairs  appears  to  raise  the  shear  stress  thus  generating  a 
tougher  network.  However,  after  the  maximum  stress  is  reached  both  curves  become  nearly 
overlapping. 

The  deformation  in  shear  mode  is  illustrated  in  Figure  29  using  six  snapshots  taken  at 
ascending  strain  values.  At  s  =  0.7,  stripes  are  stretched  to  sigmoidal-like  shape.  As  the  network 
is  stretched  further  at  s  =  1.3,  2.0  and  2.6  they  get  straighter  and  thinner  showing  that  strands  are 
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extended  at  or  beyond  their  limits  before  breaking.  Finally,  at  e  =  3.3,  close  observation  reveals 
that  the  location  of  interfacial  failure  is  near  the  middle  of  the  lower  half  of  the  Figure  29  (e). 

Discussion 

We  found  consistently  in  both  tensile  and  shear  modes  of  deformation  that  the  ionic 
network  requires  a  larger  number  of  bonds  to  be  broken  before  the  ultimate  failure  occurs.  In 
order  to  understand  the  origins  of  this  behavior,  we  first  determined  the  positions  of  the  failing 
bonds.  These  bonds,  as  identified  in  Figure  30,  indicate  the  location  of  the  interfacial  failure.  For 
example,  the  clustering  of  these  broken  bonds  in  the  lower  portion  of  the  Figure  30  (a)  points  to 
the  location  of  failure  as  confirmed  by  Fig.  3  for  the  defonnation  of  the  standard  network  in 
tensile  mode.  Figure  30  (c)  shows  that,  in  the  case  of  ionic  network,  these  bonds  are  present  not 
only  in  the  location  of  the  failure  but  also  in  the  surrounding  area  by  a  larger  quantity  than  the 
control  network  (Figure  30  (a)).  This  result  is  even  clearer  in  shear  mode  (Figure  30  (b)  vs.  (d)). 

We  infer  that  the  integration  of  ionic  pairs  bring  about  a  network  that  can  distribute  the 
stress  more  unifonnly,  which  in  turn  results  in  improved  mechanical  properties.  In  order  to  gain 
insight  into  how  the  stress  is  distributed  more  uniformly  we  looked  at  the  growth  of  voids  that 
form  with  tensile  deformation.  It  is  an  important  phenomenon  because  the  strands  surrounding 
voids  are  where  the  stresses  concentrate.  The  snapshots  of  voids  are  shown  in  Figure  31  at 
varying  strains.  At  the  yield  point  ( e  =  0.1),  voids  are  visible  in  their  early  form.  They  appear  to 
be  more  randomly  distributed  initially  in  the  standard  than  the  ionic  network,  while  at  strains  e  > 
1.0  they  are  equally  randomly  distributed.  The  variation  in  their  distribution  is  not  expected  to 
affect  the  mechanical  properties  significantly,  because  the  difference  in  randomness  is  seen  only 
until  bond  breaking  starts.  However,  it  is  notable  to  observe  that  the  sizes  of  voids  are  smaller  in 
the  ionic  network  compared  to  the  standard  network  at  any  strain.  The  smaller  void  sizes  in  the 
ionic  network  shows  that  the  load  can  be  more  uniformly  distributed  throughout  the  network, 
hence  overall  mechanical  properties  are  improved,  in  line  with  our  analysis  of  distribution  of 
broken  bonds. 

We  had  observed  earlier  that  the  shear  stress  is  enhanced  up  to  the  point  of  the  interfacial 
failure  (Figurd  28).  After  this  point,  stresses  became  indistinguishable  between  the  two  types  of 
networks.  Unlike  this  observation,  when  a  lower  deformation  rate  (by  an  order  of  magnitude)  is 
applied  (see  Figure  7),  the  ionic  network  continues  to  bear  a  higher  load  even  after  the  point  of 
the  maximum  stress.  Knowing  that  the  two  experiments  differ  only  by  the  deformation  rate  and 
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that  the  networks  differ  only  by  the  presence  of  ionic  pairs,  this  result  indicates  that  it  is  a 
consequence  of  the  rate  dependence  of  ionic  pairs.  A  plausible  explanation  is  that  the  ionic  pairs 
that  deform  as  a  result  of  shearing  are  able  to  re-form  only  when  the  system  is  deformed  with  a 
low  strain  rate  so  that  the  strands  Figure  32  have  enough  time  to  relax.  We  infer  also  that  this 
time  dependency  refers  to  the  presence  of  a  self-healing  mechanism  that  restores  the  strength  of 
the  network. 

Concluding  Remarks 

We  presented  results  for  a  coarse-grained  model  of  cross-linked  networks  with  and 
without  ionic  bonds  incorporated.  The  MD  simulations  of  the  deformation  of  the  networks  in 
tensile  and  shear  mode  showed  that  the  model  has  captured  the  fundamental  mechanical  features 
such  as  the  yield  behavior  and  the  strain  hardening. 

One  is  bound  by  the  length  scale  that  can  be  practically  modeled,  as  an  inherent  limitation 
of  molecular  level  simulations.  This  constraint  reflected  itself  for  the  current  system  in  two  ways. 
First,  how  the  stress  is  transferred  from  the  solid  layers  to  the  network  should  be  a  function  of  the 
distance  between  the  upper  and  lower  layers.  The  interfacial  failures  that  occur  in  the  vicinity  of 
the  layers,  rather  than  in  the  middle  section  in  all  four  cases  may  be  a  consequence  of  the  stress 
build  up  in  this  region.  Secondly,  the  network  can  be  only  as  heterogeneous  as  the  system  size.  In 
fact,  the  main  reason  why  the  system  fails  only  after  very  high  strains  (c  >2)  is  that  the  crack 
initiation  can  occur  only  after  strands  are  extended  significantly.  In  a  real  material  there  would  be 
a  heterogeneous  cross-linking  (i.e.  broad  distribution  of  degree  of  freedom  of  strands  to  stretch). 
When  the  deformation  is  applied,  the  strands  with  a  low  degree  of  freedom  to  stretch  will  be 
stressed  the  most  and  break  first.  These  regions  where  such  constrained  strands  concentrate  will 
act  as  defects  that  initiate  cracks  before  the  rest  could  collectively  stretch.  From  this  point  of 
view,  the  current  model  appears  to  represent  an  ideal  network  compared  to  a  real  one.  Another 
interpretation  is  that  the  model  is  capturing  the  deformation  behavior  at  or  near  the  crack  tip  as 
suggested  previously  67-69 . 

The  primary  outcome  of  this  work  is  the  finding  that  the  integration  of  ionic  pairs 
elevated  the  level  of  stress  that  the  network  can  withstand  over  a  wide  range  of  strain  values  in 
both  tensile  and  shear  mode  of  deformation.  This  implies  that  ionic  pairs  create  a  stronger  and 
tougher  network  compared  to  the  control  network.  We  identified  two  plausible  mechanisms 
through  which  the  ionic  pairs  became  effective.  They  were  able  to  reform  after  being  displaced 
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due  to  deformation  and  they  enabled  a  more  unifonn  distribution  of  the  stress  throughout  the 
network. 

In  future  work,  we  plan  to  investigate  the  effects  of  the  cross-link  density,  content  of 
ionic  pairs  and  the  interplay  between  them. 


Figure  25.  Radial  distribution  functions  of  pairs  of  (a)  chain-chain  a n d  xlinker-xl inker 
and  (b)  ion-ion  and  layer-xlinker  type  beads. 
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( s )  behavior 


'X 

upon  tensile  deformation  at  strain  rate  of  0.6X  10'  1/x.  Standard  network  (line)  and 


ionic  network 


(circles). 


Figure  27.  Views  of  the  network  during  deformation  in  tensile  mode  at  s  =  0.0,  0.5,  1.0, 


1.5  and  2.0. 
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upon  shear  deformation  at  strain  rate  of  0.6X  10'  1/x.  Standard  network  (line)  and  ionic  network 
(circles). 


Figure  29.  Views  of  the  network  through  x-z  plane  during  shear  deformation  at  strain 
values  of  (a)  0.0,  (b)  0.7,  (c)  1.3,  (d)  2.0,  (e)  2.6,  and  (f)  3.3.  The  beads  are  colored  in  stripes  to 
help  visualize  the  shearing. 
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Figure  30.  Views  (through  x-z  plane)  of  pairs  of  beads  forming  bonds  that  eventually 
brake  after  applying  deformation.  Upper  half  shows  views  for  the  control  network  after 
deformation  in  (a)  tensile  and  (b)  shear  mode.  Similarly,  lower  half  shows  views  for  the  ionic 
network  after  defonnation  in  (c)  tensile  and  (d)  shear  mode.  The  lines  highlight  approximate 
locations  of  ultimate  interfacial  failure. 
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Figure  3 1 .  Snapshots  of  the  voids  as  they  develop  during  tensile  deformation  in  the 
standard  (upper  panel)  and  the  ionic  (lower  panel)  networks.  Strains  at  which  snapshots  were 
generated  are  (from  left  to  right):  0.1,  0.15,  0.25,  0.5,  1.0  and  1.5. 


Figure  32.  Stress  (cr)  versus  strain  (;;)  behavior  upon  shear  defonnation  at  strain  rate  of 
0.6X10’4  1/t.  Standard  network  (line)  and  ionic  network  (circles). 
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2.4  Mechanisms  of  Gecko  Adhesion 

Gecko’s  exceptional  ability  in  climbing  up  smooth  vertical  surfaces  has  intrigued  many 
scientists  leading  to  a  detailed  investigation  of  the  mechanisms  of  gecko’s  adhesion.  Microscopy 
has  shown  that  a  gecko's  foot  has  a  hierarchy  of  structures:  nearly  five  hundred  thousand  30-130 
pm  long  setae  in  the  foot  branching  into  hundreds  of  projections  terminating  in  0.2 -0.5  pm 
spatula-shaped  structures  (see  Figure  33,  left).  While  the  strong  adhesion  ability  is  attributed  to 
this  hierarchical  structure  and  the  rapidly  releasing  ability  to  the  anisotropy  in  the  spatula,  there 
is  no  fundamental  understanding  of  the  basic  adhesion-deadhesion  mechanism.  The  knowledge 
of  this  mechanism  would  have  important  implications  in  designing  responsive  adhesive 
biomimetic  systems.  We  plan  modeling  the  unique  structure  of  the  spatula  as  a  networked  system 
in  a  coarse-grained  representation  of  the  repeating  unit  (Figure  33,  right).  We  plan  to  compare 
the  simulation  results  to  the  experimental  measurements  of  the  directional  adhesion  of  model 
systems  and  to  bring  insight  into  the  fundamentals  of  the  process. 

Thus  far,  we  used  a  two-dimensional  coarse-grained  approach  to  model  a  cross-linked 
parallelepiped  body  with  a  50°  tilt.  This  shape  is  intended  to  mimic  the  type  of  anisotropy  that 
appears  to  be  the  underlying  mechanism  in  gecko’s  adhesion-release  ability  to  surfaces. 
Computer  simulations  of  detachment  experiments  at  varying  detachment  directions  were  carried 
out.  Early  results  have  shown  a  qualitative  agreement  with  experimental  measurements  of 
directional  adhesion  for  a  similar  system  by  Guduru  at  Brown  University. 


Figure  33.  (Left)  Hierarchical  structure  of  gecko’s  foot.  (Right)  A  model  of  tilted  pillar 
in  coarse-grained  representation. 
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TIMELINE  OF  TECHICAL  PROGRESS 
December-2005 

In  the  past  month  I  have  reviewed  the  current  literature  on  the  newly  assigned  project  of 
developing  ice-phobic  materials.  I  have  worked  on  introducing  the  idea  of  using  the  anti-freeze 
proteins  for  developing  ice-phobic  materials.  It  has  been  recently  discovered  that  antifreeze 
proteins,  found  in  organisms  that  can  survive  in  extremely  cold  environments  (e.g.  Antarctic 
fish),  disrupt  ice  growth  in  a  non-colligative  manner.  These  proteins  may  directly  be  utilized  as 
ice  inhibitors  or  the  knowledge  of  how  they  function  can  be  used  for  designing  novel  coating 
systems.  I  am  planning  on  modeling  the  coating  surface-  anti-freeze  protein  interactions.  This 
computational  effort  may  potentially  initiate  an  innovative  path  towards  developing  ice-phobic 
materials. 

I  have  also  researched  computer  hardware  and  software  that  will  be  required  for  the 
modeling  and  simulation  effort.  Particularly,  I  have  configured  and  collected  quotes  for  a 
workstation  and  a  Beowulf  cluster.  As  configuring  Beowulf  clusters  require  a  great  deal  of 
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knowledge  in  concepts  such  as  parallel  computing,  message  passing  interface,  network 
interconnect  etc.  I  have  spent  considerable  amount  of  time  understanding  these  concepts  and 
configuring  the  optimum  hardware  for  the  task. 

January-2006 

In  the  past  month  I  have  worked  on  solidifying  the  project  outlines  for  the  ice-phobic 
materials  project  by  reviewing  the  literature  on  the  modeling  of  ice  surface.  For  that  purpose  I 
have  analyzed  various  types  of  crystal  forms  of  ice.  I  built  the  unit  cell  for  hexagonal  ice  (Ih),  so 
that  as  soon  as  the  necessary  AMBER  modeling  program  is  installed  on  the  workstation  I  can 
replicate  the  unit  cell  to  build  the  simulation  box.  This  ice  structure  will  be  essential  for  the 
future  simulation  work,  where  the  interactions  between  a  candidate  ice-phobic  material  and  the 
ice  surface  will  be  investigated. 

I  also  worked  with  Capt  Kris  Hardy  on  a  LDF  proposal.  This  proposal  is  based  on  the 
modeling  and  simulation  of  the  ice-phobic  materials  project.  Specifically,  the  adhesive 
interactions  between  ice  and  other  surfaces  will  be  investigated.  These  interactions  between 
certain  chemical  groups  are  proposed  to  provide  a  guideline  in  designing  ice-phobic  materials  in 
the  laboratory.  Additionally,  modeling  the  binding  of  ice  to  an  antifreeze  protein  will  be 
performed.  This  will  provide  an  insight  into  how  the  ice  fonnation  problem  is  solved  in  nature 
which  in  turn  will  provide  clues  as  to  how  this  problem  may  be  solved  via  synthetic  methods. 
Therefore,  how  they  function  can  be  used  to  design  biomimetic  coating  systems  with  inherent 
anti-icing  capability. 

February-2006 

In  the  past  month  I  dedicated  a  considerable  part  of  my  time  to  setting  up  my  Linux 
workstation.  I  worked  with  a  computer  support  specialist  for  installation  of  Linux  operating 
system  as  well  as  the  following  tools:  Intel  Math  Kernel  Libraries,  FFTW  libraries,  MPICH 
parallel  processing  libraries,  Intel  C/C++  compilers,  Intel  Fortran  compiler,  AMBER  molecular 
dynamics  simulator,  LAMMPS  molecular  modeling  software,  VMD  visualization  software, 
SPDBV  modeling  tool  and  XMGRACE  scientific  graphics  software.  I  have  tested  all  of  this 
software.  Particularly,  I  modeled  a  water  box  of  -600  water  molecules  using  Tip4p  and 
Tip4p/Ice  water  potentials  and  ran  a  test  simulation  at  300  K.  The  purpose  of  this  modeling  and 
simulation  effort  was  two  folds.  First,  I  checked  the  AMBER  for  correctly  running  and 
benchmarked  its  perfonnance  on  the  workstation.  Secondly,  I  wanted  to  evaluate  the  qualities  of 
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the  water  potentials  for  reproducing  certain  properties  such  as  density.  Tip4p/Ice  water  potential 
is  particularly  important  as  it  has  been  shown  to  reproduce  crystallization  into  ice  lh  phase  at  a 
temperature  very  close  to  the  experimental  observation.  Therefore,  I  expect  this  potential  to  be  at 
utmost  importance  for  the  ice-phobic  materials  research. 

I  also  improved  the  ice  lh  crystal  I  created  previously.  In  addition  to  the  positions  of  the 
oxygen  atoms  that  has  to  satisfy  the  experimentally  determined  unit  cell  for  ice  lh,  the  positions 
of  the  hydrogen  atoms  is  also  important.  Namely,  there  is  randomness  in  the  hydrogen  bonding. 
This  randomness  has  to  satisfy  the  “ice  rules”  detennined  by  Bernal  and  Fowler.  Additionally, 
the  unit  cell  has  to  have  zero  net  dipole  moment  and  minimal  net  quadrupole  moment.  I 
implemented  the  results  of  a  work  by  Hayward  and  Reimers  so  that  all  these  requirements  are 
met  without  having  to  go  through  time  consuming  Monte  Carlo  simulation  to  create  the  unit 
cells.  This  implementation  is  realized  by  a  Fortran  code  I  developed  using  GNU  Fortran77.  The 
resulting  unit  cells  can  directly  be  used  for  any  future  modeling  studies  that  involve  ice  lh. 

A  tertiary  effort  that  I  engaged  in  the  past  month  was  to  acquire  the  best 
performance/price  for  a  Beowulf  cluster  that  we  would  like  to  purchase  for  the  modeling  work  in 
AFRL/MLBT.  As  a  result  of  this  effort  we  agreed  on  a  deal  with  an  established  vendor  that 
provides  us  with  the  most  economical  package  with  equivalent  product  specifications  and  with 
on  site  maintenance  that  other  vendors  would  not  provide.  I  am  confident  that  the  modeling 
effort  at  MLBT  will  benefit  substantially  from  the  acquisition  of  this  Beowulf  cluster. 
March-2006 

In  the  past  month,  I  have  modeled  a  polyurethane  (the  main  constituent  in  topcoat  of  an 
aircraft  paint)  slab  in  order  to  create  the  surface  on  which  the  ice  will  adhere  and  the  adhesive 
interactions  of  will  be  studied.  The  molecule  used  in  the  modeling  is  a  dimer  of  a  product  of  the 
reaction  of  m-Phenylene  diisocyanate  with  ethylene  glycol.  This  molecule  has  been  chosen  for 
its  simplicity  in  order  to  evaluate  the  success  of  the  force  field  and  the  simulation  protocol. 
While  AMBER  molecular  dynamics  program  package  has  shown  that  it  is  capable  of  modeling 
and  simulating  the  amorphous  polyurethane  system  at  the  realistic  density  successfully,  it  lacks 
some  features  that  enable  conveniently  modeling  of  a  two  phase  system  such  as  ice  and 
polyurethane.  Therefore,  my  current  effort  is  to  incorporate  LAMMPS  molecular  modeling 
package  for  modeling  the  interface.  I  developed  programming  codes  in  Fortran  90  for  the 
following  purposes.  First  codes  that  make  AMBER  coordinate  and  topology  files  portable  to 
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LAMMPS  environment  and  vice  versa.  Then  programs  that  analyze  output  from  LAMMPS 
simulations  and  calculate  averages  and  standard  deviations  for  several  quantities  and  plot  the 
results  in  graphs  automatically. 

A  second  Linux  workstation  to  be  used  in  simulations  has  recently  become  available.  As 
I  have  done  for  the  first  workstation  I  have  been  working  with  a  computer  support  specialist  for 
installation  of  Linux  operating  system  as  well  as  the  following  tools:  Intel  Math  Kernel  Libraries, 
FFTW  libraries,  MPICH  parallel  processing  libraries,  Intel  C/C++  compilers,  Intel  Fortran 
compiler,  AMBER  molecular  dynamics  simulator,  LAMMPS  molecular  modeling  software, 
VMD  visualization  software,  SPDBV  modeling  tool  and  XMGRACE  scientific  graphics 
software.  I  also  worked  networking  two  workstations  so  that  the  files  will  be  transferred  easily 
back  and  forth.  The  acquisition  of  this  second  computer  will  reduce  the  overall  simulation  time. 
April-2006 

I  have  modeled  an  isocyanurate  molecule  (Figure  a)  and  an  oligomer  of  polyester  polyol 
(Figure  b).  The  crosslinking  reactions  of  the  two  molecules  produce  crosslinked  polyurethane; 
the  type  that  is  commonly  used  in  aircraft  topcoats.  I  have  then  developed  FORTRAN  code  that 
creates  a  simulation  box  of  binary  blends  on  a  cubical  lattice.  Using  this  initial  configuration 
created  by  this  program  I  have  performed  equilibration  simulations  of  isocyanurate  and 
polyester.  Using  this  equilibrated  structure  I  have  built  another  simulation  box  of  the  topcoat 
(isocyanurate  and  polyester)  and  the  ice  crystal.  This  initial  configuration  was  further 
equilibrated  and  the  ice-topcoat  interface  was  obtained  (Figure  c).  Next,  I  will  begin  working  on 
calculating  the  energy  of  adhesion. 

For  the  purpose  of  making  crosslinked  polyurethane  thus  modeling  of  a  more  realistic 
coating  system,  I  have  also  been  working  on  the  modeling  the  crosslinking  reaction.  As  reactions 
are  not  supported  by  simulation  software,  it  was  necessary  for  me  to  develop  a  programming 
code  that  performs  this  task.  I  have  completed  writing  a  code  that  can  react  given  two  atoms. 
Now,  I  will  be  working  on  extending  the  code  so  that  it  can  perfonn  a  search  in  the  simulation 
box,  detennine  the  pairs  of  the  closest  reacting  groups  and  then  perfonn  the  crosslinking 
reaction. 
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May-2006 

Adhesion  of  Ice  on  Coating  Surfaces 

I  have  completed  developing  the  crosslinking  code  “xlinker”.  The  resulting  over  1300 
line  long  code  will  likely  be  an  indispensable  tool  for  creating  cross! inked  coating  systems 
towards  the  modeling  and  simulation  efforts  in  the  paints  and  coatings  group.  The  main  function 
of  the  code  is  to  determine  the  candidate  reacting  pairs  and  carry  out  the  reactions.  In  doing  this 
it  sorts  the  reacting  pairs  with  respect  to  their  closeness  thus  closer  pairs  are  given  priority  in  the 
order  of  reaction.  From  the  modeling  point  of  view,  the  reaction  of  a  pair  of  reacting  groups 
mean  that  a  set  of  changes  occur  on  several  atoms  in  the  reacting  molecules.  These  modifications 
include  change  of  atom  types,  bond  types,  angle  types  and  dihedral  angle  types.  Additionally, 
some  existing  bonds,  angles  and  dihedrals  removed  and  new  ones  created. 

An  important  feature  of  the  code  is  that  it  perfonns  crosslinking  reaction  of  molecules  in 
the  primary  simulation  box  with  their  images.  Therefore,  one  can  build  crosslinked  molecules 
with  infinite  length  in  any  of  the  three  dimensions.  The  modeled  system  therefore  becomes  very 
similar  to  the  real  systems.  Figure  (a)  shows  a  two  dimensional  projection  of  a  crosslinked 
molecule  that  is  infinite  in  x  and  y  directions. 

The  crosslinked  system  needs  to  be  characterized  in  order  to  know  the  micro-scale 
architecture  of  the  crosslinked  molecules.  The  code  performs  a  complete  analysis  and  reports  a 
wide  set  of  information  including  density  of  crosslinking  reactions,  percentage  of  reactions, 
number  of  newly  formed  chains  and  their  size,  the  members  of  newly  formed  chains,  which 
reactive  group  in  which  molecules  is  linked  with  which  other  groups  etc. 

Currently  I  am  running  simulations  of  crosslinked  polyurethane-ice  interface.  As  these 
simulations  finish  I  will  extract  the  adhesion  energies  between  polyurethane  and  ice. 

Tunable  Adhesion 

We  have  recently  initiated  a  project  on  the  adhesion  strength  and  mechanisms  of  single 
polymer  chains  with  polar  end  groups  on  various  surfaces  within  the  context  of  tunable  adhesion. 
The  computer  modeling  of  the  system  by  MD  simulation  enables  to  obtain  the  force- 
displacement  curve  which  then  helps  us  to  understand  the  non-covalent  interactions  between 
surfaces  and  candidate  molecules  for  tunable  adhesion.  This  modeling  effort  will  be  carried  out 
along  with  scanning  probe  microscope  (SPM)  experiments.  While  the  SPM  experiments  helps  to 
validate  the  simulation  results,  the  computer  modeling  enables  to  perfonn  studies  of  several 
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molecules  with  varying  architectures  and  end  groups  that  are  not  available  commercially  and 
very  difficult  to  synthesize  in  the  lab.  Therefore,  the  computer  simulation  and  the  SPM 
experiments  will  synergistic  ally  provide  some  insight  into  the  molecular  interactions  leading  to 
adhesion. 

I  have  built  a  pair  of  gold  crystal  surfaces.  One  of  them  is  the  substrate  and  the  other  one 
is  the  SPM  probe.  I  adapted  a  new  potential  function  along  with  its  force  field  (Dreiding  FF)  in 
order  to  correctly  predict  gold-organic  molecule  interactions.  I  have  run  a  simulation  in 
LAMMPS  where  the  gold  substrate  was  kept  stationary  whereas  the  probe  was  displaced  at  a 
constant  speed.  The  polyethylene  chain  with  thiol  groups  at  the  end  was  sandwiched  between  the 
gold  surfaces.  This  configuration  (Figure  (b))  mimics  SPM  experiments  so  that  the  results  of  the 
simulation  can  be  compared  with  the  experiments.  My  current  effort  is  to  calculate  the  force  on 
the  golden  probe  that  is  imposed  by  the  polymer  chain  as  the  probe  moves  upwards. 

June-2006 

I  have  continued  working  on  the  tunable  adhesion  project.  In  addition  to  the  single  chain 
configuration,  I  have  created  a  self  assembled  monolayer  (SAM)  configuration.  Accordingly,  I 
modeled  a  SAM  made  up  of  1 ,4-Benzenedimethanethiol  (BDMT)  and  then  placed  it  between  a 
couple  of  gold  substrates  (see  Fig.  1).  Based  on  experimental  and  theoretical  evidence  I  placed 
thiol  groups  on  a  (V3XV3)-R30°  triangular  lattice  on  a  surface  of  Au  (111).  The  preliminary 
equilibration  simulations  of  the  modeled  system  showed  that  the  BDMT  molecules  oriented  as  to 
from  a  herring-bone  pattern  as  shown  in  Fig.  2  (inline  with  experimental  evidence).  Next  effort 
will  be  to  perform  the  simulations  of  AFM  experiment  on  the  system  and  calculate  force- 
displacement  curve. 

I  have  authored  computer  modeling  part  of  an  AFOSR  proposal  based  on  the  tunable 
adhesion  project.  Basically,  we  proposed  to  reveal  adhesion  mechanisms  by  combining  data  from 
the  atomic  force  microscope  (AFM)  and  the  molecular  dynamics  simulations. 

I  have  attended  a  three  day  workshop  (2006  AFMC  Deicing  Workshop)  held  in  WPAFB. 
I  prepared  a  report  based  on  the  workshop  to  be  submitted  to  my  government  sponsor. 

July-2006 

I  have  continued  working  on  the  tunable  adhesion  project.  I  performed  multiple  computer 
simulations  of  Atomic  Force  Microscope  (AFM)  experiment  for  1 ,4-Benzenedimethanethiol 
(BDMT)  in  the  single  molecule  and  the  self  assembled  monolayer  (SAM)  configurations.  I 
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obtained  the  force-distance  curves  which  is  the  main  objective  of  this  work.  Currently,  I  am 
waiting  for  input  from  the  actual  AFM  experiments  to  compare  the  simulation  results  to  those 
obtained  from  the  experiments. 

I  have  worked  on  incorporating  the  visualization  of  the  computer  simulations  to  the 
evolution  of  the  force-distance  curve  in  the  real  time.  Namely,  my  objective  has  been  to  display 
how  the  force  imposed  on  the  AFM  tip  by  the  BDMT  molecules  progresses  as  a  result  of 
approaching  the  tip  towards  the  BDMT  molecule/s.  For  this  purpose  I  taught  myself  a  scripting 
language  (TCL)  that  can  do  the  job  efficiently  and  easily.  The  resulting  TCL  code  inputs  the  data 
from  the  visualization  software  (VMD)  and  stream  it  to  the  to  the  graphics  plotting  tool 
(XMGRACE)  in  the  real  time.  The  code  then  converts  the  collection  of  resulting  images  into  a 
movie  format  (mpg)  that  can  be  displayed  on  any  type  of  computer.  The  outcome  of  this  effort 
enables  presenting  the  simulation  results  clearly  by  combining  the  qualitative  results  to  the 
quantitative  ones  in  the  real  time.  Snapshots  taken  from  the  created  movies  for  the  single 
molecule  and  SAM  configurations  are  shown  in  Figures  1  and  2,  respectively. 

August-2006 

I  have  continued  working  on  the  tunable  adhesion  project.  Previously,  I  had  written  a 
“tel”  script  that  calculates  the  three  angles  (9,  \|/,  (j>)  that  detennine  the  conformation  of  BDMT 
molecules  and  the  distribution  functions  for  each  of  these  three  angles.  In  an  attempt  to  carry  out 
similar  analyses  for  other  self  assembling  molecules,  I  modeled  biphenyldithiol  (BPDT), 
benzylmercaptan  (BM),  4-nitrobenzenethiol  (NBT),  thiophenol  (TP),  4-methylbenzenethiol  (TT) 
and  4-methylbenzylthiol  (XT).  I  am  currently  working  on  obtaining  the  equilibrium 
configurations  of  the  self  assembled  monolayers  (SAMs)  of  these  molecules. 

Another  point  of  interest  for  this  study  is  the  determining  the  positions  (fee,  hep  or 
bridge)  of  sulfur  atoms  at  the  bottom  of  self  assembling  molecules  on  the  Au  (111)  surface.  For 
this  purpose  I  wrote  a  “tel”  script  that  determines  the  position  of  each  of  270  sulfur  atoms  on  the 
gold  surface  and  tracks  it  throughout  the  course  of  the  molecular  dynamics  simulation.  Having 
that  information,  with  the  FORTRAN  code  I  wrote,  one  can  determine  the  correlation  function 
for  each  possible  site  (fee,  hep  or  bridge).  This  function  shows  on  average  how  much  time  a  self 
assembling  molecule  spends  on  each  of  the  site. 
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September-2006 

I  have  continued  working  on  the  tunable  adhesion  project.  Previously,  I  had  written  a 
“tel”  script  that  calculates  the  three  angles  (9,  \| /,  (j))  that  detennine  the  conformation  of  BDMT 
molecules  and  the  distribution  functions  for  each  of  these  three  angles.  In  an  attempt  to  carry  out 
similar  analyses  for  other  self  assembling  molecules,  I  modeled  benzylmercaptan  (BM),  4- 
nitrobenzenethiol  (NBT),  thiophenol  (TP),  4-methylbenzenethiol  (TT)  and  4-methylbenzylthiol 
(XT)  (see  Fig  1).  I  am  currently  working  on  obtaining  the  equilibrium  configurations  of  the  self 
assembled  monolayers  (SAMs)  of  these  molecules. 

Another  point  of  interest  for  this  study  is  determining  the  positions  (fee,  hep  or  bridge)  of 
sulfur  atoms  at  the  bottom  of  self  assembling  molecules  on  the  Au  (111)  surface.  For  this 
purpose  I  wrote  a  “tel”  script  that  determines  the  position  of  each  of  270  sulfur  atoms  on  the  gold 
surface  and  tracks  them  throughout  the  course  of  the  molecular  dynamics  simulation.  Having  that 
information,  with  the  FORTRAN  code  I  wrote,  one  can  determine  the  correlation  function  for 
each  possible  site  (fee,  hep  or  bridge).  This  function  shows  on  average  how  much  time  a  self 
assembling  molecule  spends  on  each  of  the  site  (see  Fig  2). 

I  also  wrote  a  “tel”  script  to  calculate  the  distance  between  the  sulfur  atom  at  the  bottom 
of  self  assembling  molecules  and  its  3  nearest  gold  atoms.  I  compared  this  data  to  the  results 
from  Density  Functional  Theory  methods  and  observed  a  very  good  correspondence. 
October-2006 

I  have  been  working  on  the  “adaptable  adhesion”  project  with  an  emphasis  on  modeling 
the  self-assembled  monolayers  (SAMs).  In  addition  to  previously  modeled  SAM  of  1,4- 
Benzenedimethanethiol,  I  modeled  and  obtained  the  equilibrium  configurations  of  SAMs  for 
benzylmercaptan,  biphenyldithiol,  thiophenol,  4-methylbenzenethiol  and  4-methylbenzylthiol  on 
Au  (111)  surface.  I  have  performed  various  analyses  on  the  SAMs  in  order  to  characterize  and 
understand  their  conformational  and  dynamical  characteristics.  Specifically,  I  completed  writing 
codes  and  scripts  to  obtain  quantities  including: 

•  distribution  functions  of  three  angles  ( 6,  y/,  (f>)  that  determine  the  confonnation  of 
self-assembling  molecules 

•  positions  (fee,  hep  or  bridge)  of  sulfur  atoms  of  SAMs  on  the  Au  surface 

•  time  dependent  correlation  functions  for  each  possible  site  (fee,  hep  or  bridge)  (see 
Figure  1) 
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.  distance  between  the  sulfur  atom  of  self  assembling  molecules  and  its  3  nearest  gold 
atoms 

.  mean  square  displacement  (msd)  for  self-assembling  molecules  (see  Figure  2) 

Based  on  the  outcome  of  these  analyses,  I  am  currently  working  on  writing  a  paper  to  be 
submitted  to  a  scientific  journal. 

I  also  attended  a  3 -day  workshop  on  scientific  software  development  at  Ohio 
Supercomputer  Center  (Columbus)  in  October  24  -  26,  2006. 

November-2006 

I  have  been  working  on  writing  an  article  based  on  my  modeling  results  of  the  self- 
assembled  monolayers  (SAMs).  I  completed  “model”  and  “simulation  protocol”,  sections  while 
“results”  and  “discussion”  sections  are  near  complete. 

I  have  also  been  working  on  configuring  the  new  Beowulf  cluster  that  I  now  have  access 
to.  I  am  working  with  the  cluster  vendor  and  the  computer  support  in  installing  necessary 
software  and  resolving  various  configuration  issues. 

An  additional  effort  was  on  the  call  for  Statement  of  Need  on  “Scientific  Understanding 
of  Non-chromated  Corrosion  Inhibitors  Function”  by  Strategic  environmental  Research  and 
Development  Program  (SERDP).  I  have  done  preliminary  literature  search  on  the  matter  and 
proposed  some  ideas  to  approach  the  problem. 

December-2006 

I  have  been  working  on  writing  an  article  based  on  my  modeling  results  of  the  self- 
assembled  monolayers  (SAMs).  I  now  completed  the  introduction,  results  and  discussion 
sections  in  addition  to  the  model  and  simulation  protocol  sections.  I  am  planning  to  complete  the 
paper  within  this  coming  month. 

I  worked  on  composing  two  LDF  (Laboratory  Director’s  Fund)  proposals  with  Dr.  Alex 
Khramov. 

First  proposal  is  titled:  “Controllable  nano-structured  interface  based  on  self-assembled 
monolayer  and  gold  nano-particles”.  We  proposed  to  create  a  hierarchical  structure  with 
controllable  interfacial  properties  by  sandwiching  a  SAM  between  a  gold  surface  and  a  layer  of 
gold  nano-particles.  Our  approach  is  based  on  versatile  and  well-developed  chemistry  of  thiol- 
functionalized  SAMs  on  a  gold  surface  and  thiol-stabilized  gold  nano-particles.  In  addition  to  the 
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experimental  approach,  we  plan  to  develop  atomistic  models  of  an  ordered  and  a  disordered 
SAM  sandwiched  between  Au  (1 1 1)  layers. 

Our  second  proposal  is  titled:  “ Surface  modification  and  patterning  by  grafting 
molecular  loops  of  alkylene-bridged  bis-silanes” .  The  goal  of  this  proposal  is  to  achieve 
effective  surface  modification  of  a  metal  substrate  by  anchoring  alkylene-bridged  bis-silanes  to 
the  surface  thus  creating  interfacial  structure  that  consists  of  molecular  loops  extended  out  of  the 
surface.  Moreover,  we  will  perform  molecular  dynamics  simulations  in  order  to  elucidate  the 
thennodynamic  as  well  as  the  structural  mechanisms  involved  in  this  surface  modification 
process. 

January-2007 

I  have  completed  the  article  based  on  my  modeling  results  of  conformational  and 
dynamical  behaviors  of  the  self-assembled  monolayers  (SAMs).  It  is  currently  being  internally 
reviewed.  I  plan  to  submit  it  to  Langmuir. 

Our  proposal  with  Dr.  Alex  Khramov  entitled:  “Controllable  nano-structured  interface 
based  on  self-assembled  monolayer  and  gold  nano-particles  ”  qualified  for  the  second  step  of  the 
selection  process.  Thus,  I  worked  on  composing  an  extended  proposal  that  was  re-submitted  to 
the  LDF  (Laboratory  Director’s  Fund)  committee. 

In  addition  to  my  previous  simulations  for  ordered  SAMs,  I  modeled  a  disordered  layer 
of  benzenedimethanethiol  (BDMT)  on  a  Au(l  1 1)  surface.  Subsequently,  I  performed  simulations 
of  AFM  experiment  on  the  disordered  surface.  The  analyses  for  this  simulation  results  will  reveal 
the  role  of  ordering  in  adhesive  function  of  self  assembling  molecules  between  Au  layers. 
February-2007 

I  initiated  a  new  project  on  modeling  of  self-healing  polymers.  The  first  stage  of  the 
project  involves  using  a  course-grained  approach  in  modeling  the  cross-linked  polymer  system.  I 
developed  FORTRAN  codes  that  fonn  the  initial  structure  of  reactants  which  will  fonn  the  cross- 
linked  system.  Additionally,  I  am  working  on  developing  a  scheme  that  dynamically  cross  links 
the  reactants  by  using  home  developed  x-linker  code  and  communicating  with  LAMMPS 
molecular  simulator. 

I  started  preparing  my  presentation  to  be  presented  at  ACS  National  Meeting  on  March 
26,  2007. 
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I  attended  a  Python  Programming  workshop  on  21  February,  2007  at  Ohio  Super 
Computer  Center. 

March-2007 

I  continued  the  project  on  self  assembled  monolayers  as  a  part  of  more  general  effort  on 
tunable  adhesion  project.  I  perfonned  substantially  longer  10  ns  long  simulations  for  six  aromatic 
thiol  molecules.  I  recalculated  some  quantities  such  as  mean  squared  displacement  and  site- 
correlation  functions  with  the  larger  set  of  data  (Fig.  1).  I  incorporated  this  data  into  the  paper 
that  is  to  be  sent  to  Langmuir  after  our  internal  review. 

I  continued  the  work  on  a  course-grained  approach  in  modeling  the  crosslinked  system 
for  self-healing  project.  I  improved  the  FORTRAN  code  that  form  the  initial  structure  of 
reactants  which  subsequently  forms  the  crosslinked  system  and  the  scheme  that  dynamically 
cross  links  the  reactants  by  using  the  x-linker  computer  code  and,  by  communicating  with 
LAMMPS  molecular  simulator.  The  next  step  is  to  perfonn  uniaxial  stretching  simulations  and 
to  calculate  the  stress-strain  curves. 

I  attended  a  5  day  long  2007  American  Chemical  Society  National  Meeting  in  Chicago, 
where  I  presented  a  talk  entitled  “Molecular  Dynamics  Simulation  Study  of  Conformational  and 
Dynamic  Properties  of  Self- Assembled  Thiol  Monolayers  on  Au”. 

April-2007 

I  continued  the  work  on  the  course-grained  approach  in  modeling  the  crosslinked  system 
for  self-healing  project.  The  FORTRAN  code  that  forms  the  initial  structure  of  reactants  and 
subsequently  forms  the  crosslinked  system  was  improved.  I  updated  the  computational  scheme 
that  dynamically  cross  links  the  reactants  by  using  the  computer  code  x-linker  and,  by 
communicating  with  LAMMPS  molecular  simulator.  Subsequently,  I  perfonned  simulations 
where  three  modes  of  deformations  (i.e.  tensile,  shear  and  compression)  were  applied  to  the 
polymeric  control  network.  The  stress  response  to  the  various  modes  of  defonnation  was  probed, 
which  completes  the  first  stage  of  the  project.  The  next  step  will  be  the  incorporation  of  ionic 
cross  linking  points  to  the  control  network  (i.e.  creating  ionic  network).  The  comparisons  of  the 
stress-strain  characteristics  for  the  control  network  and  the  ionic  network  will  be  utilized  to 
quantitively  assess  the  self-healing  potential  of  the  ionic  network. 
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May-2007 

I  have  worked  on  a  talk  was  presented  at  CERMACS  2007  May  20-23,  2007  entitled: 
“Exploring  conformations  and  dynamics  of  self-assembled  thiol  monolayers  on  Au  by  molecular 
dynamics  simulations.” 

I  have  worked  on  an  abstract  for  an  invited  talk  to  be  presented  at  SAMPE  Fall  Technical 
Conference,  October  29-November  1,  2007  entitled:  “Adding  Self  Healing  Character  to 
Polymeric  Networks  via  Ionic  Bonds:  A  Coarse-Grained  Molecular  Dynamics  Simulation 
Study.”  This  talk  will  be  based  on  our  work  where  we  adapted  a  coarse-grained  representation  in 
modeling  of  the  network  structure  which  was  achieved  by  dynamically  crosslinking  the 
reactants.  It  involves  our  results  for  the  stress-strain  curves  after  imposing  various  forms  of 
deformation  and  measuring  their  stress  responses. 

I  also  worked  on  finalizing  the  journal  article  that  I  submitted  to  Langmuir  after  the 
internal  review. 

I  have  been  doing  technical  supervising  for  a  summer  hire  on  the  self-assembled 
monolayer  research. 

June-2007 

I  have  been  working  towards  a  paper  to  be  submitted  to  SAMPE  for  the  Fall  Technical 
Conference.  This  paper  entitled  “Adding  Self  Healing  Character  to  Polymeric  Networks  via 
Ionic  Bonds:  A  Coarse-Grained  Molecular  Dynamics  Simulation  Study.”  is  based  on  the  work 
where  a  coarse-grained  representation  was  adapted  in  modeling  of  the  network  structure  which 
was  achieved  by  dynamically  crosslinking  the  reactants.  For  this  purpose,  I  have  been 
performing  simulations  of  various  forms  of  deformation  and  measuring  their  stress  responses  in 
order  to  obtain  the  stress-strain  curves.  I  have  also  been  doing  calculations  on  various  properties 
based  on  the  simulation  results. 

I  have  been  doing  technical  supervising  Sheena  Neely  on  the  self-assembled  monolayer 
research. 

I  worked  with  computer  support  in  configuring  one  of  the  Linux  workstations.  We 
installed  the  operating  system,  compilers  and  various  modeling  and  simulation  software  and 
tools. 
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July-2007 

We  are  initiating  a  collaborative  work  with  Brian  Hinderliter  at  North  Dakota  State 
University  on  modeling  and  simulation  of  coating  systems.  As  an  initial  step,  I  started  reviewing 
the  literature  on  appropriate  atomistic  force  fields  and  coarse-graining  a  force  field  for  polymer 
networks. 

I  have  been  researching  software  that  performs  quantum  chemistry  calculations  and  that 
visualizes  the  models  and  the  outputs  from  the  calculation.  We  acquired  Gaussion  03  quantum 
chemistry  software.  I  also  evaluated  “Molekel”  and  “gOpenMol”  visualization  tools.  The 
objective  is  to  add  quantum  level  (ab  initio )  modeling  capability  to  our  modeling  resources. 

I  have  been  supervising  the  summer  researcher  Sheena  Neely  as  a  technical  point  of 
contact  on  the  self-assembled  monolayer  (SAM)  research.  I  helped  her  in  building  the  models, 
running  the  simulations,  performing  the  analyses  on  linear  alkanethiol  SAMs.  I  also  advised  her 
during  preparation  of  a  technical  paper  and  a  presentation. 

August-2007 

I  have  been  working  on  coarse-grained  modeling  of  highly  cross-linked  polymeric 
networks  towards  the  self-healing  polymeric  networks  project.  The  previous  system  size  was 
extended  from  4,000  beads  to  100,000.  This  increase  aims  to  obtain  a  more  representative  and 
accurate  model  as  well  as  less  scattered  data  for  properties  such  as  stress-strain  curves  and 
toughness.  The  cross-linking  procedure,  which  I  previously  designed,  results  in  a  very  high 
cross-linking  density  of  about  98%.  In  addition,  the  degree  of  ionic  pairs  was  increased  by  50%. 
I  am  currently  running  equilibrium  simulations  of  the  system  that  will  be  followed  by 
deformation  simulations. 

Another  focus  of  my  efforts  has  been  obtaining  the  force-distance  curves  for  self 
assembled  monolayers  (SAMs).  I  intend  to  model  SAMs  of  benzenebiphenyldithiol  (BPDT)  at 
three  different  possible  phases:  ordered  standing  up,  ordered  lying  down  and  disordered.  In  my 
previous  work  I  had  modeled  the  ordered  standing  up  phase  correctly.  Therefore,  I  have  been 
working  on  modeling  the  ordered  lying  down  phase.  As  experimental  data  on  this  particular 
phase  is  missing,  its  modeling  has  been  intricate.  I  have  recently  managed  to  resolve  the  relevant 
issues  and  obtained  the  model.  I  am  currently  working  to  model  the  disordered  phase. 
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September-2007 

I  have  continued  working  on  coarse-grained  modeling  of  highly  cross-linked  polymeric 
networks  towards  the  self-healing  polymeric  networks  project.  After  building  of  a  system  with 
100,000  beads  and  obtaining  equilibrium  structure,  I  have  been  working  on  optimizing  the 
simulation  parameters  for  deformation  of  the  system  in  different  modes.  The  test  simulations 
have  shown  that  the  new  set  of  parameters  produce  stress-strain  curves  in  tensile  and  shear 
modes  with  desired  consistency  and  accuracy  (see  Figure).  In  addition,  another  set  of  systems 
with  varying  degree  of  ionic  cross-linking  points  have  been  created.  The  objective  is  to 
understand  the  effect  of  ionic  cross-linking  density  on  the  mechanical  and  self-healing 
properties.  I  am  currently  working  on  a  script  that  will  run  a  series  of  simulations  automatically 
on  a  supercomputer. 

October-2007 

I  have  worked  on  completing  the  simulations  for  the  self-healing  polymeric  networks 
project.  Stress-strain  curves  were  obtained  after  exposing  the  standard  model  and  ionic  model  to 
three  different  modes  of  deformations  at  higher  strain  rates.  The  results  showed  the  expected 
enhancements  in  mechanical  properties,  particularly  in  impact  resistance,  were  observed.  The 
simulations  for  slow  deformation  rate  are  being  run. 

I  worked  on  preparing  a  talk  in  Computational  Materials  Science  session  at  SAMPE  Fall 
Technical  Conference  that  I  delivered  on  Oct  3 1st. 

I  am  currently  working  on  an  invited  talk  at  the  Department  of  Chemical  &  Materials 
Engineering  in  University  of  Cincinnati  on  Nov  8,  2007. 

November-2007 

I  have  been  working  on  creating  self-assembled  monolayers  of  biphenyldithiol  (BPDT)  in 
lying-down  as  well  as  standing-up  configurations.  In  addition  to  the  ordered  monolayers  I  have 
created  a  disordered  thin  film  of  the  same  compound.  I  have  perfonned  computational  AFM 
experiments  to  obtain  the  force-distance  curves  on  the  three  different  phases.  The  objective  is  to 
compare  the  force-distance  behavior  for  the  different  surfaces. 

I  have  presented  an  invited  talk  at  the  Department  of  Chemical  &  Materials  Engineering 
in  University  of  Cincinnati  on  Nov  8,  2007. 
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I  have  authored  a  short  paper  entitled  "Course-grained  molecular  dynamics  simulations  of 
self  healing  polymeric  networks",  for  a  presentation  at  6th  International  Conference  on 
Mechanics  of  Time-Dependent  Materials. 

December-2007 

I  have  continued  working  on  studying  self-assembled  monolayers  of  biphenyldithiol 
(BPDT)  in  lying-down  as  well  as  standing-up  configurations.  The  computational  AFM 
experiments  have  shown  that  the  structure  of  monolayer  has  a  significant  impact  on  the  adhesive 
character  of  the  surface. 

I  have  started  writing  a  full  paper  to  be  submitted  to  the  Mechanics  of  Time  Dependent 
Materials  Journal.  This  paper  will  be  based  on  the  coarse-grained  molecular  dynamics 
simulations  of  self-healing  polymeric  networks. 

I  have  also  pursued  an  exploratory  research  to  understand  the  effect  of  directional 
anisotropy  on  adhesion.  I  used  a  coarse-grained  approach  to  model  a  cross-linked  body  with  a 
unique  shape  that  is  expected  to  provide  anisotropy. 

January-2008 

I  have  continued  an  exploratory  research  to  understand  the  effect  of  directional 
anisotropy  on  adhesion  that  is  seen  in  gecko  foot.  I  used  a  two-dimensional  coarse-grained 
approach  to  model  a  cross-linked  parallelepiped  body  with  a  50°  tilt.  This  shape  is  intended  to 
mimic  the  type  of  anisotropy  that  appears  to  be  the  underlying  mechanism  in  gecko’s  adhesion- 
release  ability  to  surfaces.  Computer  simulations  of  detachment  experiments  at  varying 
detachment  directions  were  carried  out.  Early  results  have  shown  a  qualitative  agreement  with 
experimental  measurements  of  directional  adhesion  for  a  similar  system  by  Guduru  at  Brown 
University. 

I  have  completed  writing  a  full  paper  to  be  submitted  to  the  Mechanics  of  Time 
Dependent  Materials  Journal.  This  paper  is  based  on  the  coarse-grained  molecular  dynamics 
simulations  of  self-healing  polymeric  networks  and  it  addresses  the  effect  of  addition  of  ionic 
pairs  into  the  network  on  time  dependent  mechanical  properties. 

February-2008 

I  have  been  working  on  organizing  the  results  for  the  study  of  computational  AFM 
experiments  on  self-assembled  monolayers  of  biphenyldithiol  (BPDT).  The  results  of  force- 
distance  measurements  in  standing-up,  lying-down,  disordered  configurations  are  being 
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compiled  for  a  manuscript.  An  additional  case  when  the  AFM  tip  was  coated  with  residual 
dithiols  prior  to  the  approach  and  withdrawal  stages  in  the  AFM  simulation  was  studied.  I  have 
started  writing  the  manuscript  for  a  possible  submission  to  Langmuir. 

I  have  modeled  another  cross-linked  network  system  for  the  self-healing  polymeric 
networks  project  in  order  to  validate  our  earlier  results  of  stress-strain  curves  after  applying 
different  modes  of  deformation. 

I  provided  input  for  the  TRB  review  of  “Adaptive  &  Active  Composite  &  Hybrid 
Materials  Program”.  I  summarized  the  computational  modeling  effort  and  the  current  status  in 
RXBT  in  regards  to  the  adaptable  adhesion  and  self-healing  materials  projects.  In  addition,  I 
compiled  the  input  provided  from  other  researchers  working  on  other  aspects  of  the  program  to 
be  submitted  in  the  write-up  for  the  review. 
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